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ABSTRACT 


Classically, the T-matrix method is a procedure to exactly compute the multiple scattering of an 
incident wave from a “cloud” of objects, given knowledge of the free-held scattering properties 
of a single object for an arbitrary incident wave. For acoustic waves. Profs. Baker and Scan- 
drett have extended the T-matrix method to the case in which the radiation sources are also the 
scatterers, that is, to the case of an array of active transducers. This thesis is the hrst successful 
practical demonstration of the T-matrix method applied to an active sonar array for which a 
hnite-element model was employed to compute the scattering properties of a single transducer. 
For validation, a T-matrix model of a linear array of piezoelectric spherical thin-shell transduc¬ 
ers was modeled, for which analytical approximate values of the T-matrix element values are 
known. Subsequently, a T-matrix model of a linear array of piezoelectric class V flextensional 
“ring-shell” transducers was modeled. Beam patterns of the linear array models computed with 
the T-matrix method are compared with those of an array of point sources, demonstrating that 
the T-matrix method produces more realistic beam patterns, especially for end fire arrays. 
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CHAPTER 1: 

Introduction 


1.1 Motivation 

In many sonar applications, it is desirable to use a high power and directional sonar transducer. 
Their use reduces reverberation and can propagate further distances. Also, they provide more 
aeeurate estimate of bearing. These transducers are used in arrays in order to maximize the 
directionality and steer optimum acoustic beams. Figure 1.1 shows two examples of volumetric 
sonar array applications. 



(a) (b) 

Figure 1.1: (a) A bi-planar array eontaining of 24 Sparton of Canada 600 Hz elass V flex- 
tensional transducers [1] (b) Bicycle Array composed two of the 610 Hz class V flextensional 
transducers [2] 

The interaetion between array elements can affect array performance. Especially, this effect 
becomes severe for high power and directional array elements because they are commonly 
operated at a resonance frequency and are spaced closely, typically no greater than A/2 [8]. 
Aeeurate performanee modeling of such an array requires proper accounting of the acoustic 
interaction between array elements. 

Historically, the Pritchard’s approximation [9] has been employed to evaluate inter-element ar¬ 
ray interactions. It has been employed to model a bi-planar array by the Naval Underwater 
Systems Center (NUSC, New London, CT; now the Naval Undersea Warfare Center, NUWC, 
Newport, RI) to predict the accurate beam pattern. Beam pattern predictions of computer mod¬ 
els were then compared against results of deep water testing. Although the eomputer models 
resulted in mueh better beam patterns with regard to the non-interacting point sources model, it 
showed a significant deviation at and near end fire steering angles [8]. 
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Another example of evaluating array interaction comes from a fundamental acoustic course at 
the Naval Postgraduate School (NPS). It is a three-element array with omnidirectional drivers 
that was investigated for several different configurations in an anechoic chamber. The investiga¬ 
tion uncovered substantial deviation between theoretical and experimental results. In particular, 
near end fire, source-level amplitude data does not match theory, and the product theorem [10] 
cannot explain this discrepancy Figure 1.2. 


Polar Plot for 
847 Hz 



EXP. 

THEO. 


0 


Figure 1.2: The far-field beam pattern for a three-element array with omnidireetional drivers [3] 

In the early 1990s, a "T-Matrix Method" was developed by NPS researchers Prof. S. R. Baker 
and Prof. C. L. Scandrett [6,11] in order to properly model randomly and densely packed arrays 
of active sonar transducers. 

To build an array model using the T-Matrix method, one must obtain the free-field radiation 
properties for a given electrical drive condition and the scattering properties (T-Matrix coeffi¬ 
cients) for a single-transducer, from such as a finite-element model (FEM). 

This thesis focuses on building a three-dimensional (3D) FEM of a volumetric piezoelectric 
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transducer for T-matrix calculation in COMSOL and properly developing a model of a closely 
packed linear array by using the T-Matrix Method. 

1.1.1 Related Works 

Although the T-matrix method was first developed in the 1990s, it has not been applied to a 
realistic volumetric array model because, until now, the capabilities and user-friendliness of the 
finite-element computer codes available for performing the single-transducer radiation and scat¬ 
tering computations required for the T-matrix method were insufficient. Thus, the method was 
not of practical use to a non-specialist engineer. The attempt of prior theses to compute T-matrix 
elements using ATILA finite-element code resulted in excessive discrepancies in comparison to 
analytical results [12,13]. 

Recently, that situation has changed. Using the COMSOL multiphysics finite-element computer 
code, finite-element model computations for the acoustic scattering T-matrix coefficients of a 
thin-wall elastic spherical shell now show excellent agreement with those computed using full 
3D elasticity theory [14]. 

Therefore, the T-matrix method is applicable and can be used to build an exact model of densely 
and randomly-packed active sonar transducer array. Furthermore, this method can achieve a 
more accurate analysis of array performance that better defines the acoustic characteristics of 
the array. 

1.2 Purpose 

The purpose of the research described in this thesis is to demonstrate for the first time the 
practical application of the T-matrix method to the numerical performance modeling of an ac¬ 
tive sonar array, that is, to combine the results of finite-element model computations of the 
single-transducer free-field radiation and scattering characteristics with an analytic model for 
the free-field propagation between transducers. Specifically, the T-matrix method will be ap¬ 
plied to model the source level and far-field beam pattern for several example unsteered linear 
active sonar arrays comprised of either 1) radially-polarized piezoelectric spherical thin-shell 
transducers or 2) type V flextensional "ring-shell" transducers, model 34A0610, manufactured 
by Sparton of Canada [7,15]. The former was chosen as an academic example, as approx¬ 
imate analytic radiation and scattering properties for such transducers are known; the latter 
was chosen because it represents a real high-power, low-frequency sonar projector which was 
built and tested in a 24-element, densely-packed, bi-planar, active array, and for which single- 
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element and array measured performance characteristics are available [7,16]. Although only 
unsteered linear array configurations of 1/2-wavelength spacing and "ka"~l were considered 
herein, the ultimate goal of the research program, of which the present effort is a contribu¬ 
tion, is to compare the characteristics of the as-tested bi-planar array as computed using the 
T-matrix method to both its measured characteristics and its characteristics as computed using 
the so-called "Pritchard approximation" to model array element interaction [9]. 

1.3 Organization 

The remainder of the thesis consists of six chapters. 

Chapter 2 clarifies features of COMSOL’s Acoustic Module and presents the addition theorem 
used in coordinate translation matrices for an array. The second part of chapter 2 describes the 
far-field beam pattern calculation for a linear array using the T-matrix method. 

Chapter 3 describes the characteristic features of the modeled transducers which consist of a 
radially polarized spherical transducer and a type V Flextensional Ring-Shell Model 34A0610. 

Chapter 4 is devoted to building 3D models of the radially polarized piezoelectric spherical thin- 
shell transducer and the type V flextensional ring-shell transducer. The second part of Chapter 
4 includes an evaluation of sought-after transducer parameters from the 3D FEM model. 

Chapter 5 calculates the total acoustic pressure field and produces the far-field beam pattern 
plots for equally spaced linear arrays using the T-matrix methodology. 

Chapter 6 concludes this thesis and provides suggestions for future work. 

1.4 Assumptions 

Single-frequency harmonic time dependence of the form is assumed throughout this thesis. 
The symbol is used for the imaginary unit to avoid confusion with the spherical Bessel 
function of the first kind, 
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CHAPTER 2: 
Theory 


2.1 COMSOL Finite Element Analysis 

2.1.1 COMSOL Multiphysics 

The COMSOL Multiphysics numerical package is a finite element code that includes simula¬ 
tion software capable of modeling many different physical phenomena employing a variety of 
different physics-based modules. The modules involve a solution method of Partial Differential 
Equations (PDE) and predefined multiphysics interfaces. Predefined multiphysics interfaces 
simplify the model by using a built-in coupled physics node. They undertake the analysis for 
different types of studies (e.g., modal, eigenfrequency, frequency response) [17]. 


2.1.2 The Acoustic-Piezoelectric Interaction Interface 

The Acoustic-Piezoelectric Interaction Interface is one of the physical modes of COMSOL 
Multiphysics within the Acoustics Module. This module can solve the far-field or near-field ra¬ 
diation problem by incorporating acoustic-structure interactions of piezoelectric materials [18]. 
The boundary conditions include the pressure load from fluid to solid and a structural accelera¬ 
tion between the solid and fluid. The Acoustic-Piezoelectric Interaction Interface combines the 
piezoelectric device, pressure acoustics, frequency domain, solid mechanics, and electrostatic 
interfaces [18]. 


2.1.3 Acoustics Interface Theory in COMSOL Multiphysics 

The acoustic wave equation is expressed in COMSOL by the following equation [18]: 

Bulk 

Modules 
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Variable 

Description 

Unit 

P 

Fluid Density 

kg/m^ 

c 

Speed of Sound 

m/s 

p 

Pressure 

Pa 

fld 

Dipole Source 

N/m^ 

Qm 

Monopole Source 

l/s2 


In the case of plane waves, the pressure varies with time as 

p{x,t) = p{x)e'^\ (2.2) 

The assumption of the frequency domain or time-harmonic dependence for the acoustic wave 
equation results in the inhomogeneous Helmholtz equation [18]: 

= (23) 

Boundary conditions are required to solve the inhomogeneous Helmholtz equation. COMSOL 
provides several different types of boundary conditions, depending upon the analysis performed. 
The prescribed boundaries between the fluid-solid boundary are added after nodes are selected 
and fixed. Three built-in boundary conditions are available [18] on the boundaries between the 
fluid and solid: 

1. free 

2. sound hard boundary(wall) 

3. acoustic structure boundary 

A free boundary condition assumes the medium reacts without any constraints acting on the 
boundary. 

A sound hard boundary is added for a rigid boundary. With this boundary condition, the normal 
component of differential velocity is zero [18]. This boundary condition is also named the 
Neumann condition: 

n- (^(Vp-q^)^ =0 (2.4) 

An acoustic structure boundary couples the fluid load and structural acceleration at the fluid-solid 
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boundary. It cannot be added to other boundaries [18]. 


¥p = -np 

= n • Uft 


(2.5) 


Variable Description 

F Force 

p Pressure 

n Unit normal,outward pointing from solid domain 

an Normal acceleration for the acoustic fluid 

utt Second derivative of the Structural Displacement with respect to time 


2.1.4 Piezoelectric Device Theory in COMSOL Multiphysics 

Piezoelectric materials are solids for which the mechanical strain and the electric field in the 
solid influence each other. Changes in strain or electric field cause responses of the polarization 
and stress within the solid. For example, an applied potential difference induces a deformation 
of the piezoelectric device, and/or the deformation of the piezoelectric material induces an 
electrical potential difference within the piezoelectric device. The piezoelectric relations for 
linear systems having small strains are expressed by the following equations in COMSOL [18]. 

Strain-Charge : 

S = SET + d^E ( 2 . 6 ) 

D = rfT + £rE (2.7) 

This relation can be rewritten into stress-charge form as 
Stress-Charge: 

T = C£S - (2.8) 

D = eS + £5E (2.9) 

Piezoelectric ceramics are “transverse isotropic” materials [19]. In material science, orthogonal 
axes, designated “1,” “2,” and “3” are defined by convention for a given material type. The 
polarization is conventionally designated as the “3” direction, and COMSOL uses the same 
convention for piezoelectric devices. 
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The piezoelectric material properties are expressed in stiffness, compliance, coupling, and di¬ 
electric matrices with respect to their local coordinate system. The default local coordinate 
system corresponds to the Cartesian coordinate system [17]. 


Variable 

Description 

T 

Stress 

S 

Strain 

D 

Electrical displacement 

se 

Elastic compliance constant 

CE 

Elastic stiffness constant 


Dielectric constant 

e 

Piezoelectric stress constant 

d 

Piezoelectric strain constant 


2.2 Piezoelectric Spherical Thin-Shell Transducer 

The spherical thin shell transducer is a particular transducer in which the spherically symmetric 
mode (lumped constant mode) can be developed from an equation of motion rather than a wave 
equation (Figure 2.1) [4,20]. Therefore, fundamental acoustic features (e.g, mechanical Q, 
radiation load) of a thin shell transducer can be accurately represented by a simple closed form 
expression. The equivalent circuit [4] represents a model of a thin shell transducer (Figure 2.2). 
Hence, it is particularly convenient for illustrating equivalent circuits and design principles [20]. 

A radially polarized spherical thin shell transducer of mean radius “a”, and thickness “t” is 
modeled. The radial direction coincides with the material ”3“ direction. Electrodes are located 
on the inner and outer surfaces of the sphere to form equipotential surfaces. Hence, transverse 
fields E\ and E^, which are perpendicular to the radial direction, are neglected. All of the tensile 
stresses are equal to zero except the circumferential tensile stress 7) [4,19]. The relationship 
between strain and stress forces in the “1” and “2” directions are obtained from symmetry 
The relationship between strain and stress forces in the 1 and 2 directions are obtained from 
symmetry 

Ti=T2 = Tc 

S\=S2= Sc- 

Compliance coefficients are combined to form a circumferential compliance coefficient ex- 
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O V Q 



Figure 2.1: Radial cross section of a radially polarized piezoelectric spherical thin shell 
transducer 


pressed as 

= 2 (^ 11 +‘^ 12 ). 


(2.10) 

The piezoelectric equations are then 


Sc = 2 s^Tc + d 3 iE 3 

(2.11) 

D 3 = 2^3 1 Tc + & 33 E 3 . 

(2.12) 


The electrical displacement equation can be expressed by substituting in equation (2.11) 


D3 = 


-^Sc + £ 33 ( 1 — kp^) £'3. 


(2.13) 


Where 


- r^T 




(2.14) 


The circumferential strain equals the ratio of the radial displacement to mean radius of the 
sphere 

Sc = -. (2.15) 

a 


The radial displacement can be obtained from the equation of motion, which is derived from 
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Newton’s Second Law 


where, 


The equation of motion becomes 


Fr + pwta56^-7r = 0 
at‘^ 

Fr = F8e = 2TcWt8e. 


c 

pwta8d^:^ = —2TcWt8d. 

oF 


(2.16) 


The equation of motion can be expressed using differentials by substituting Equations (2.11) 
and (2.15) into Equation (2.16) 


1 e <^31 ^ 


aps^ 


dt^ a^ps^ 

The free resonance frequency of the undriven sphere (E=0) can be written from 2.17 

2 _ 1 

The spherical shell displacement for harmonic excitation is then expressed as 

•,2T7^Ja)t ^^..r.^ 2 ] 


(2.17) 


(2.18) 




adi,\ (O^Eoe^^^ ad^iCO^E^ 


co^-co^ 


(0^-(0^ 


(2.19) 


After using the relationship between symmetrical strain and spherical shell displacement, Sc = ^ja, 
this equation can be written in a more convenient form for computing the current, I = Ako^D , 


i)3=£3(^| + 43(l-*:/) 


SO 


/ — ATta^icoEj, ^ + ^33 (l “ ] • 

The electrical admittance can be obtained from equation (2.21) 

I 1 ATta^iO) (^31 


( 2 . 20 ) 


( 2 . 21 ) 


m 


V z 


(©2-0)2) 


+4(1-^/) 


( 2 . 22 ) 
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The equivalent circuit for a spherical ring, Figure 2.2, can be obtained for the electrical admit¬ 
tance. 



Figure 2.2: Equivalent circuit for a spherical transducer [4] 


for which the electrical element values are: 


M = Ana^tp 

^ We|3(l-fe^) 


rE = ^ 
4m 


7 _ ^ cjkg_ 

S = 4na^ 


(2.23) 


Element 

Description 

Go 

Electrical loss conductance 

Co 

Open-circuit capacitance 

N 

Electromechanical turns ratio 


Short-circuit compliance 

M 

Lumped mass 

Rm 

Mechanical Resistance 

Zr 

Free-field radiation impedance 


2.3 Type V Flextensional Ring-Shell Transducer 

Ring-shell flextensional transducers are desirable because of their ability to operate at low fre¬ 
quencies while generating high acoustic power. They also have higher operating depth than 
uncompensated flextensionals of similar size and weight [2]. Furthermore, they provide wide 
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bandwidth and low mechanical quality factor, Qm, which are very important characteristics for 
underwater applications. 

A composite piezoelectric ceramic and steel ring is used as a motor element. The ring is sand¬ 
wiched between two convex thin spherical shell caps Figure 2.3. The small motion of the ring 
causes a relatively large volume displacement of the shell. The ring-shell transducer is depth 
compensated and operates at low frequencies. 

The desired resonance frequency can be obtained by adjusting the piezo-ring diameter. In order 
to compensate for the depth pressure, a water bladder is used. The bladder is opened to the 
ocean and adjusts the internal pressure according to the depth pressure. [15] 



Figure 2.3: Cross Sectional View of Type V Flextensional Ring-Shell Transducer 


2.4 T-Matrix Method 

2.4.1 Introduction 

Application of the T-Matrix method requires computation of the single-transducer scattering 
properties for an arbitrary incident wave. Since the incident pressure distribution is not known a 
priori, it, along with the single-transducer free-held radiation, is parametrized as a series expan¬ 
sion with to-be-determined coefficients using the free-held eigenfunctions of the linear wave 
equation in spherical coordinates. Outgoing waves (free-held radiated and scattered) are rep¬ 
resented using outgoing traveling wave solutions (spherical Hankel functions), while incident 
waves are represented using standing wave solutions (spherical Bessel functions). Unnormal¬ 
ized spherical harmonic functions are used in all cases as an angular basis set, as shown in 
equations (2.24) and (2.25) 

oo n 

p‘(r,e,^)=Y, E (2.24) 

n m=—n 
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(2.25) 


oo n 

p"(r.e,<.) = £ £ P’„„hif>(kr)a:{e,4,), 

n m=—n 

where 



P”{cosd)e’”‘l’. 

Variable 

Description 


The spherical coordinates 

pi 

^ nm 

Incident pressure coefficient 

^nm 

Scattered coefficient from the obstacle 

jn{kr) 

order spherical Bessel function 

h^n\kr) 

order spherical Hankel function of the second kind 


Unnormalized spherical harmonic function 

p“(cose) 

Associated Legendre function of the first kind 


To simplify further development, the following matrix notation will be employed. The spherical 
harmonic component coefficients for any quantity belonging to a single radiator are arranged in 
a column matrix, denoted using bold capital letter, e.g., by “flattening out” the double sequence 
over n and m, as follows 


/ ps \ 

^00 ^ 


f ps ^ 

ps 

M-1 


Pl 

ps 

MO 

= 

ps 

M 

^11 



V / 


V / 


The one dimensional running index of the column matrix equals 


1 +n +n + m 


In principle, the length of such matrices may be infinite; however, in a practical computation, 
the number of spherical harmonics kept is truncated to a maximum value up to a certain degree, 
typically the second degree (quadrupolar, n=2) [6]. If the maximum spherical harmonic degree 
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2 

kept is Umax, the total number of spherical harmonic components kept is TV = (umax + 1) > e.g., 
9 for Hffiax — 2. 


2.4.2 Unnormalized Versus Normalized Spherical Harmonics 

The unnormalized spherical harmonic function (0,0) is related to the normalized spherical 
harmonic function F™ (0,0) [21]; 

77(0,0) =C„„Q“(0,0) =C„„P“(cos0)exp(/m0) (2.26) 


where 


Cn 


(-I)" 


2n +1 {n — m)\ 
Alt (n + m)! 



(2.27) 


Any function can be decomposed into a superposition of either normalized or unnormalized 
spherical harmonic functions [21] 

oo n 

/(0,0)=E I a™ir(0,0) (2.28) 

n=0yn=—n 

where the constant is given by 

K In 

anm = jI [77(0,0)]7(0,0)r/a = C™/™ (2.29) 

0 0 


so 

oo n 

/(e,'«= E E (2.30) 

n=0m=—n 

The orthonormality property of the associated Legendre equation is used for the limit 0 to 
0-dependent and 0 to 2k 0-dependent. 

The Inm is given by [14] 

71 271 

I„m = j I [a-(0,<^)]V(0,0)^^^- (2.31) 

0 0 
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2.4.3 T-Matrix Calculation 


The multiple scattering of acoustic waves by transducers in an array, to all orders, is rigor¬ 
ously and self-consistently accounted for by enforcing that the scattered component amplitudes 
at each transducer are related to the free-field incident component amplitudes from all other 
transducers by a transition, or T-, matrix 


ps ^ _ 


(2.32) 


Employing the matrix representation introduced in Section 2.4.1, this relationship can be ex¬ 
pressed as 


' P(' 


' Tn . 


' Pi' 

P^N 


. Tnn 


ph 


(2.33) 


r\ 

where = (nmax +1) • 

The T-matrix elements are computed through a succession of scattering computations, as fol¬ 
lows. An incident wave is imposed with unit amplitude in exactly one member ( component) of 
the incident eigenfunction expansion, and the amplitudes of all members (components) of the 
scattered eigenfunction expansion found. This determines the member of exactly one column 
of the T-matrix. The T-matrix elements of that column equal the component pressure ampli¬ 
tudes. For example, to determine the first column of the T-matrix, elements denoted Rnm below, 
an incident wave with = 1 and all other = 0 is imposed, and the resulting left hand 
side column matrix of scattered wave component coefficients determined by applying or¬ 
thogonality of the unnormalized spherical harmonic functions (see Section 2.4.2 for details). If 
(a, 6, (l>) is the scattered pressure amplitude at radius a in the fluid, then 


where 


I 

(ka) 


K 2k 

Inm = j I 
0 0 


(2.34) 


(2.35) 
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ps 

^00 


ps 

^ nm 


Rqq 


R„ 


where 


0 


(2.36) 


Rnm — Pnm ~ 


I 

(ka) 


(2.37) 


This follows from setting 


oo n oo n 

E E o?’(Mn:'(0.'C) = E E ci/,„£i“(e,<.) 

n m=—n n=0M=—n 


(2.38) 


It should be noted that by its spherical symmetry property, the T-matrix for the spherical thin- 
shell transducer must be diagonal. Hence, the scattering cannot mix components. This is not 
exactly the case for the ring-shell transducer. By its axial symmetry property, scattering by the 
ring-shell transducer can mix “n” values, but not “m” values. 


2.4.4 The Radiation Pressure Coefficient Calculation 

In the absence of scattering, the free-held radiated pressure is expressed with the spherical 
Hankel function and unnormalized spherical harmonic function 

oo n 

P* (r. e, <.) = £ £ P«„hP (kr) n" (0, <■) . (2.39) 

n=0ni=—n 


The radiated pressure coefficients can be computed using the FEM. After comparing the spher¬ 
ical harmonic function expansion (2.28), the radiation pressure amplitude can be expressed in 
terms of 4^ as 


/^2 T 

pR _ '^nm^nm 

^ nwi 


{ka) 


(2.40) 


The right-hand side variable seems the same as the scattered pressure amplitude, but the surface 
integral (4^) in this case will be different as scattering has been neglected. 
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2.4.5 The Addition Theorem 


In the case of an array, a single common coordinate system is inconvenient in applying boundary 
conditions. The spherical addition theorem must be used in this situation. It provides for a 
representation of a pressure field relative to one coordinate system in terms of another. 

Consider two active spherical sources of radius “a” and separated by distance “J”, as shown 
in Figure 2.4. The outgoing pressure of sphere two in mode {n, m) can be determined as a 
superposition of free-field incident pressures on sphere one in modes {s, u): 

K“>?{kr 2 )p:;(ih,h) = YJ.p‘sJ'Hkc,)p;(euh), ( 2 . 41 ) 

5 U 

where 

(Km) = IG 12 ] (K%) (2.42) 

and G \2 can be written from [5] 

n-\-s 

E a(s,t,n,u,m)hy>{kri2)Sir“(eil,itil2) , ( 2 . 43 ) 

t = \n — s\ 
t > \m — u\ 


where 


y, 



Figure 2.4: Two coordinate frames with arbitrary relative positions and orientations from [5] 
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a{s,t,n,u,m) = ” 


(25'+ \ )(2t + \){s — u)\{t — m + u)\{n + m)\[{s-\-t + n)/'2\\ 1 
[{n + t-s) /2]! [(n + 5-?)/2]! [(5 + r-n)/2]!(5 + ? + n+l)! J 


1+max 

E (-1) 

w=wmin 


S + t — n 

{s + t — n) /2 + w 


( S H — t 

(5 + W - y/2 - M + W 


t Tl -\ -5 

{t-\-n — s)/2-\-m — u-\-w 


with 


a 

b 


a\ 

b\{a-b)\ 


2w mm = max {n —s — t,s —n — t — 2m + 2u,t — s — n-\- 2u) 


2wniax = min {s-\-s — n,n-\-t — s — 2m + 2m, n-\-s — t-\- 2m) . 


Variable Description 

ri 2 the inter-element spacing distance,“rf.” 

012 the polar angle of origin “1” with respect to origin “2” 

< 1>12 the azimuthal angle of origin “1” with respect to origin “2” 

For a line array along the polar axis, which is all that will be considered herein, 0i2 equals either 
0 or Tt, and ^12 is “0”. 

2.4.6 Total Acoustic Pressure Calculation 

The multi-scattering of acoustic waves must be considered in the acoustic field calculation. The 
acoustic field is regarded as a superposition of outgoing waves from each transducer in the 
sonar array using the eigenfunctions of the wave equation in spherical coordinates as a basis set 
(e.g., spherical Hankel functions and spherical harmonics) centered at each array element. Out¬ 
going wave components account for both the free field radiation, when a transducer is electri¬ 
cally energized, and the scattered field. 
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pout _ pR^pS 


( 2 . 44 ) 


“Self consistency” is imposed, i.e., the incident pressure on one transducer is due to the outgoing 
pressure from all other transducers. This leads to a coupled set of linear equations in terms of 
the radiated and scattered pressure coefficients. For example, for two elements 


(pf) =|TiGi2] (pJ + p|) 

and 

(pf) =|T2G2,] (Pf + Pf) . 

Variable Description 

Scattering pressure amplitude for transducer 
Glq Translation formula from system q to system / 

T/ Scattering pressure amplitude for transducer 

Free-field radiation pressure amplitude for transducer 


In the matrix form (2.45) and (2.46) can be expressed as 


I -TiGi2 
T 2 G 21 I 


pS 

1 - 

pS 


TiGi2Pf 

T2G2lPf 


a = 


The “a” matrix is defined as 

0 T1G12 

T2G21 0 

After substitution into (2.47) we get the simpler matrix form: 


\l-a] 


pS 

^2 


^ ' = a 


Pf 

Pf 


or equivalently 


pS 1 =[!-«] 


pf 

pf 


This result is easily generalized to the case of “N” transducers. The result is 


(2.45) 

(2.46) 


(2.47) 


(2.48) 


(2.49) 


(2.50) 
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/ Pi’* \ 


\PS / 


(I-a) 


/ P? \ 


VPn / 


From this equation the scattering pressure amplitudes can be found. 

The summation of the scattered and radiated pressure for each individual coordinate system of 
the transducer provides the total pressure at any field point with respect to a global origin in 
spherical coordinates. 


^max n 

p(i,e,^.) = ££ £ (2-5i) 

l=l n=0 m=-n 

where r/, 0;, and (j)i are the range, polar angle, and azimuth angle, respectively, of the field point 
with respect to the local coordinate axes of transducer In all of the arrays considered herein, 
the local coordinate axes of all the transducers are assumed the same; only the origins differ in 
location. Hence, 0/ = 0 and 0/ = 0 for all 


2.4.7 Far Field Beam Pattern 

A sonar array can be represented as a collection of arbitrarily-randomly and densely-packed 
acoustic transducers. An example for a sonar array is shown in Figure 2.5, a bi-linear array 
spherical thin shell transducer elements. The horizontal far-field beam pattern of this sonar array 
is represented in Figure 2.6 [6]. The amplitude of pressure is given in terms of a directional 
factor, H{^), and the pressure amplitude along the acoustic axis of source. Pax, 

P{r,e,<j>)=Pax(r)H(e,(l>) . (2.52) 


The pressure amplitude change is proportional to 1 /r in the far field, so equation (2.52) becomes 


limp(r,0,</)) 

r—>oo 


Pax{l)H{d,(l)) 

r 


(2.53) 


where 

^ax(l) = limp(r,0,0)r 

r^oo 

is the far-field, on-axis pressure amplitude, “referred” to a range of 1 m . 
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ORIENTATION OF SPHERICAL TRANSDUCERS 



Figure 2.5: Orientation of the spherical transducers in the 16 element array 


Equation (2.51) can be simplified using the asymptotic form of the spherical Hankel function, 
which can be derived from the standard Hankel function [22] 

In the far field, the asymptotic form of the Hankel function is inversely proportional to the 
square root of the distance 


(kr) = 

r^co n+2 V Ttkr 




(2.55) 


So, the asymptotic spherical Hankel function can be written as 


r^oo 


71 


lim {kr) = J—^I—e-‘{kr-n^-^) 


2 kr V Ttkr 


After simplification it becomes 


lim {kr) = — 

r—>00 


(2.56) 
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Figure 2.6: Source level of array, ka=l, analytically scaled by far-field of a single spherical 
transducer [6] 


Applying this form to equation (2.51), and noting that “kr” above represents k-f, the total 
acoustic pressure for the far-field limit can be obtained 


Nx g-ik-{r-roi) «max n p 

/=i^-(^-^0/) «=o m=-n^ 

{1.51) 

where the direction of k equals the direction of r and roi is the vector position of the local origin 
of axes for transducer with respect to the global origin. 


In the far field, all radial vectors ir — roi) are approximately parallel. The distance to the far 
field does not vary much for each element, so the factor roi can be neglected compared to r 
in the amplitude factors (denominators). Although the path length difference can be ignored 

for amplitude, it has to be taken into account for phase. Factoring out the common phase —kr, 

—> 

the factor k ■ roi is the phase difference between the global origin and the acoustic origin of 
transducer /. If the acoustic center of transducer / is located at {xoi,yohZoi)^ 


k-roi = kxXoi + kyjoi + kzZoi 
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which can be written as 


k-foi = k(xoisin6cos(j)+yoisin6sin(l> + zoicos6) . (2.58) 

After the substitution of this equation into equation (2.57), the far-field total acoustic pressure 
equation is obtained 


Pffir^QA) = 


^-ikr Nx ^ 


kr 


Jk-roi 


/=! 


^max 

E- 

n=0 


(«+i) 


^Imn + /^5Jn”(9,(f) (2.59) 


with k • foi given by (2.58). 
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CHAPTER 3: 

Spherical and Ring-Shell Transducer Descriptions 


The radially polarized spherical thin-shell transducer and the class V flextensional ring-shell 
transducer are picked as example applications for the T-matrix calculation of a volumetric array 
transducer. 


3.1 Piezoelectric Spherical Thin-Shell Transducer 

The piezoelectric ceramic thin shell is polarized in thickness (radial) direction, and electrodes 
are installed on the inner and outer surface of the thin shell to provide radial polarization. The 
shell is made of lead zirconate titanate (PZT-8). This material is treated as a member of the 
Navy Type III. The Navy Type III is defined as “hard,” which delineates a relative difficulty 
of depolarization of the ceramic [19]. PZT-8 is preferred for high powered transducers. The 
properties of PZT-8 are obtained from the COMSOL material library. 

The thin shell is surrounded by a water. The inner domain of the thin shell is defined as air. 
The frequency of the transducer is defined with respect to “A:a=l,” which is the typical “ka” 
value of low frequency transducers [23]. Main operational and dimensional parameters are the 
following: 


Parameter 

Value 

Notation 

Frequency 

474 Hz 

/ 

Applied Electric Potential 

1 V 

y 

Shell outer Radius 

0.5 m 

a 

Shell thickness 

0.01 m 

t 

Radiation boundary dimension 

2.5 m 

R 


Table 3.1: The spherical thin-shell transducer dimensional and operational parameters 


The thin-shell transducer can be properly represented by an equivalent circuit and the radiation 
load can be calculated by again using an equivalent circuit (See for details Chapter 2.2). 


3.2 The Class V Flextensional Ring-Shell Transducer 

The class V flextensional ring-shell transducer (Sparton Model 34A0610), is shown in Fig¬ 
ure 3.1, and has been selected as a realistic model of a volumetric transducer. The ring-shell 
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transducer is versatile in the case of low-frequency and high-power applications [2], It is com¬ 
prised of a composite piezoelectric ceramic and steel ring sandwiched between two thin convex¬ 
shaped steel shell segments. The radial motion of the piezoelectric ring causes large volume 
displacement over the shell and provides high acoustic output [15]. 



Figure 3.1: The depth-compensated ring shell projector [2] 

3.2.1 Characteristics of the Transducer 

The outer piezo-ring diameter is 81 cm. The ring is made of piezoelectric material (PZT-4), 
which is a motor element of the transducer. PZT-4 is classified as a “hard” ceramic type (Navy 
Type I) and is used primarily for high-powered transducers [19]. The piezo-ring consists of 144 
plates of thickness-poled lead zirconate titanate, ceramic staves connected electrically in paral¬ 
lel. Rectangular-shaped ceramic staves are separated by 72 steel wedges to create an annular 
shape form [7]. Each steel wedge separates two bonded ceramic plates. Steel wedges are axially 
longer than ceramic staves in order to protect a ceramic surface from the stress that is caused by 
the shell edge [2]. A tension applied (25-40 MPa) fiber wrapping encloses the piezo-ring. The 
fiberglass wrapping is bonded with an epoxy resin. The fundamental operational characteristics 
of transducer are given by [1,15,24] 

The piezoelectric elements are polarized tangentially and are connected in parallel. The polar¬ 
ization direction is shown in Figure 3.2. 

In common with the flextensional transducer, the operating depth of the ring-shell projector is 
limited to a few meters without a depth compensating system. Different systems (e.g., liquid, 
compliant tubes, self-supporting, etc.) may be used to overcome the depth limited operation. 
The class V flextensional ring-shell projector is depth compensated by an internal water bladder 
that is opened to the sea through the water inlet [15]. The water bladder is installed into the 
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Parameter 

Frequency 

Maximum Source Pressure Level 

Efficiency 

Operational Depth 

Bandwidth at Resonance(-3 dB points) 


Value 
610 Hz 

213 dB re 1 juPa at 1 m (driven by 3000 Vrms) 

%90 at resonance, %65 at -6 dB 

Exceeds 400 m 

>140 Hz (470 Hz to 750 Hz) 


Table 3.2: Ring-Shell Transducer Operation Characteristic 



interior cavity of the shell segments (Figure 3.3). The air domain between the shell segments is 
compressed by the water bladder with respect to the depth pressure. The water bladder depth 
compensation system decreases the damping and supports shells against the depth pressure 
without constraining motion of the shell. 


The in-air weight without a foundation is 245 kg, while the in-water weight is 272 kg. 



Figure 3.3: Cross sectional view of ring shell projector [2] 
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CHAPTER 4: 
The 3D Models 


The numerical finite-element models of a radially polarized thin-shell spherical transducer and 
a class V flextensional ring-shell transducer were developed using the COMSOL Multiphysics 
version 4.3a finite-element code [25]. The acoustic-piezoelectric interaction, frequency domain 
interface of the acoustic module was selected using multi-physics because the fluid-structure 
coupling properties are enhanced at boundaries (See Chapter 2.1.3 for details). 

The development of models is comprised of the following steps: 

1. Selection of the COMSOL module and physics interface 

2. Adding the required type of analysis (i.e,frequency domain, eigenfrequency, frequency- 
domain modal) 

3. Creating a geometry in the selected space dimension (i.e, 3D, 2D, 2D axisymmetric) 

4. Setting material properties 

5. Assigning material geometrically 

6. Applying the physics interaction to the model 

7. Meshing 

8. Computation 

9. Postprocessing 

4.1 General Model Settings 

Development settings are the same for both models. 
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1. The built-in acoustic-piezoelectric interaction,frequency domain interface (ACPZ) was 
used. This interface couples the field boundary conditions between fluid and the piezo¬ 
electric material. The normal component of acceleration over the piezoelectric material’s 
surface generates sound within the fluid. Inversely, the piezoelectric material undergoes 
a mechanical distortion due to applied fluid acoustic pressure over its surface [18]. 

2. The free-field incident waves are referred to as ''the background pressure field” in COM- 
SOL [18]. The background pressure field is defined for each spherical harmonic compo¬ 
nent (n, m). In each computational run, exactly one incident component is applied within 
the fluid domain. The incident pressure amplitude was set to “1” so the incident standing 
wave equation (2.24) becomes 

oo n 

p‘{r,e,^) = Y. E Mkr)a'pe4). 

n m=—n 

3. The spherical wave radiation condition is applied on the outer boundary of the fluid do¬ 
main at radius 2.5 m to reduce the reflection on the boundary and to truncate the numerical 
domain. It models acoustic propagation without reflection from the boundary. 

4. The reference pressure for sound pressure level (SL) was selected as water (1/t Pa) from 
the SL settings window in COMSOL. 

5. The maximum size of mesh elements was limited with respect to the wavelength of sound 
in water. To achieve a reliable solution, the wavelength has to be resolved by a 10-12 de¬ 
grees of freedom (DOFs) corresponding to six elements in a wavelength. In other words, 
the mesh element size should be no larger than A /5 [ 18]. This criteria was applied strictly 
to each mesh sequence. 

6. The T-matrix elements were calculated for monopole, dipole, and quadrupole spherical 
harmonic components. 

7. Two studies have been examined in each model, i.e, study 1 for calculation of the 
single-element free-field radiated pressure while study 2 corresponds to determination 
of the T-matrix elements, each using the post processing feature. 
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System Requirements 

The Apple iMac MC813LL/A employed for COMSOL has the following technical specifica¬ 
tions: 

a. 2.7 Ghz Intel Core i5 processor 

b. 32 GB 1333 MHz DDR3 memory 

4.2 The Spherical Thin-Shell Transducer 

The main purpose for building a model of the radially polarized piezoelectric spherical thin-shell 
transducer is to calculate the elements of its T-matrix and the coefficients of its free-field radi¬ 
ated pressure. The validation of the COMSOL Multiphysics ACPZ interface was evaluated 
by comparing the obtained values of the free-field radiated pressure field at the surface of the 
spherical thin-shell transducer with an equivalent circuit model of a spherical transducer. 

4.2.1 Geometry 

The thickness of the model shell is 0.01 m and has a radius of 0.5 m. The piezoelectric ceramic 
is surrounded by water (Figure 4.1). The cavity of the thin-shell is air filled. 



Figure 4.1: Pictorial representation of the 3D spherical thin-shell transducer and fluid domain. 

The non-dimensional parameter ka=\ was selected to simulate the typical transducer operating 
condition and the frequency was adjusted to enforce ka=\. The model parameters can be seen 
in Table 3.1. 
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The surfaces of the spherical thin-shell were covered with electrodes. The inner surface of 
thin-shell was grounded while the outer surface was driven by a 1 V peak amplitude potential 
(Figure 4.2). 


A 1 



Figure 4.2: Electric potential of the thin-shell 


4.2.2 Material 

Three different material types were assigned to the geometry. 

Lead Zirconate Titanate(PZT-8) 

Properties of PZT-8 were assigned to the spherical shell. It is a very common piezoelectric 
material also known as a Navy Type III that is often used for power transducers. This material 
provides a maximum strength against high tensile forces. Material properties are already built 
into the COMSOL Material Library, and are shown in Table 4.1. 

The default polarization direction is defined along the “z” axis in COMSOL [18], but in the 
present case, the polarization direction is radial. In order to obtain this polarization a user- 
defined base vector system was added to the “definition” branch . 
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Parameter 


Value Unit 


Notation 


Density 

7600 

kg/m^ 

Compliance Matrix 

.E 

^11 

11.5 

lO-'^mVN 

.E 

^12 

-3.7 

lO-'^mVN 

^13 

-4.8 


^33 

13.5 


.E 

31.9 

lO-'^mVN 

.E 

^66 

30.4 

lO-'^mVN 

Coupling Matrix 

d-ii 

-97 

lO-'^C/N 

d33 

225 

lO-'^C/N 

di5 

330 

lO-'^C/N 

di5 

330 

lO-'^C/N 

Relative Permittivity 

efi/eo 

1290 


4/^ 

1000 


£0 

8.85 

lO-'^F/M 


Table 4.1: PZT -8 material properties (From the COMSOL material library) 


The spherical coordinate system was defined using the unit basis vector (i, j, k) [21] 


= Ce = cos 6 cos (^i + cos 6 sin — sin 0k 
X 2 = =-sin^i + cos^j (4.1) 

X 3 = e^ = sin 0 cos ^i + sin 0 sin + cos 0k . 

After the spherical coordinate system was added in the definition branch, it was selected as a 
coordinate system from the selection menu in the material properties window. 

The X 3 unit vector of the spherical coordinate system shows the polarization direction, and in 
this model it is referred to as the radial direction. 
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Water 

Water was used to define the fluid domain which surrounds the transducer. The material prop¬ 
erties were assigned as shown in Table 4.2. 


Parameter Value Notation 

Density 1000 kg/m^ p 

Speed of Sound 1490 m/s c 

Table 4.2: Water material properties 


Air 

Air was added to define the interior cavity of the thin-shell. The material properties are built 
into COMSOL. 


4.2.3 Meshing 

User-controlled meshing was used for the sequence type to control and optimize the numerical 
mesh. The “free mesh” algorithm is preferred over the geometry. The spherical transducer 
domain and fluid domain used a finer mesh size employing tetrahedral elements. The entire 
model consists of 421,598 elements (Figure 4.3). The average mesh quality is 0.73 where the 
quality measure is a number between 0 and 1 (Figure 4.4). An acceptable mesh quality is 
typically larger than 0.1 according to the COMSOL user’s manual [17]. 

According to the COMSOL user’s manual, the maximum element size of mesh should be 
smaller than A/5 (0.628 m in water at 474 Hz) in an acoustical model [18]. We have cho¬ 
sen the maximum element size of the mesh to be A/11.5 or 0.273 m which ensures accuracy. 
The statistical mesh quality for boundaries is represented in Table 4.3 


Boundaries 

Element Type 

Amount 

Avg. Element Quality 

Radiation 

Triangular 

2,852 

0.98 

Piezo material outer 

Triangular 

19,416 

0.99 

Piezo material inner 

Triangular 

18,504 

0.99 


Table 4.3: The meshing properties of the 3D spherical transducer model with a “finer mesh” 
size for all domains. 
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yix 


Mesh 



Figure 4.4: Mesh quality of the model 
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4.2.4 The Free-Field Radiated Acoustic Pressure Field 

The free-held radiated acoustic pressure held is plotted versus the on-axis distance from the 
thin-shell transducer surface. Figures 4.5 and 4.6 show the spatial free-held radiated acoustic 
held. Pressure magnitude, refer to one meter was found as 


P{\) = r p = 9.76 Pa.m . 
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Free-Field Radiated Acoustic Pressure (Pa) 



Figure 4.5: The free-held radiated aeoustic pressure held (magnitude) vs distance from the 
center of the spherical thin-shell transducer 


freq(l)=474 

TOTAL ACOUSTIC PRESSURE FIELD MAGNITUDE (Pa) 

A 19.519 

15 

10 
5 
0 
-5 
-10 

\ 

T 5.74X10'® 

Figure 4.6: The multi-slice plot of the free-held acoustic radiated held (magnitude) in “Pa” for 
a modeled spherical thin-shell transducer 



4.2.5 Approximate Analytical Radiation Pressure Calculation for a Thin- 
Shell Transducer 


The free-held radiated acoustic pressure can be expressed by the product of the mechanical 
radiation impedance, Z^, and the surface velocity, u, integrated over the surface of the transducer 


P = 


ZrU 

4na^ ' 


(4.2) 


The radiation impedance for the monopole sphere which has a uniform, radial velocity distri- 
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bution and uniform pressure on the surface is given as [4] 


Zy = ATta^pc 


{kaf'+ika 

\+{kaf' 


(4.3) 


and the velocity over the spherical transducer surface can be calculated by using an electrical 
equivalent circuit for the spherical transducer (Figure 2.2). In this figure of the equivalent 
circuit, the right side of the transformer contain the components of the mechanical impedance. 
Rfn is negligible compared to the stiffness and the compliance, so the mechanical impedance 
becomes approximately 


Zm = / (OM 


(OCE 


(4.4) 


with velocity 

NV 

U = 7-r , 

(Zy + Zffi) 


(4.5) 


where NV is the piezoelectric driving force. 


The lumped mass, M, the mechanical compliance, C^, and the electromechanical turns ratio, N, 
can be obtained from equation (2.23). 


The total surface acoustic radiated pressure is calculated by substitution of equation (4.3) and 
(4.5) into equation (4.2). Comparison of the results for an analytical calculation and the COMSOL 
computation are shown in Table 4.4. 


Method 

Real 

Imag 

Amp 

approximate analytical solution 

-9.6 

17.88 

20.29 

COMSOL Result 

-9.4 

17.11 

19.52 


Table 4.4: Comparison of the free-field radiated pressure at the surface of the spherical trans¬ 
ducer computed in COMSOL against those calculated from an equivalent circuit. 


The results in Table 4.4 show that the COMSOL result and the approximate analytical solution 
are in agreement. The magnitude value of the COMSOL result differs by 3.7 percent from the 
approximate analytical solution. This is an acceptable difference, since the approximate ana¬ 
lytical solution is derived from the equivalent circuit model of a spherical transducer, based on 
the assumption that its small cross-sectional dimensions are negligible when compared with its 
radius and with no internal losses [20], whereas the COMSOL model is fully three-dimensional. 
The approximate solution does provide insight as to the physical behavior of the transducer. 
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4.2.6 T-Matrix Calculation 

Nine different spherical harmonic components (n=0,l,2 m=-n..n) of an incident wave pressure 
were used in defining the background pressure field. The amplitude of a given incident wave 
pressure mode was set to “1” with the remaining mode amplitudes set to zero and scattered 
component amplitudes computed (see Section 2.4.3 for details). For example, to determine the 
first column of the T-matrix, an incident wave with Pqq = 1 and all other = 0 is imposed. 
Then resulting column matrix of scattered wave component calculated using (4.6) 


and 


I 

^nm — (4.0) 

hn (ka) 


... 


Rqq . 


1 

ps 

^ nm 


Rnm . 


0 


Evaluation of the surface integrals were performed at the surface of the spherical thin-shell 
transducer using the solution data set. One column of the T-matrix was derived from each 
computation, so nine computations were undertaken in order to complete the T-matrix. The 
T-matrix of a spherical thin-shell shell transducer is given in Table 4.5. The accuracy of the T- 


T-matrix Elements 



pi 

Mo 

pi 

M-1 

pi 

MO 

pi 

Ml 

pi 

M-2 

pi 

M-1 

pi 

Mo 

pi 

Ml 

pi 

M2 

ps 

^00 

6.91E-02 

2.60E-08 

1.18E-08 

7.36E-08 

7.51E-08 

4.51E-07 

1.94E-06 

6.99E-06 

2.81E-06 

ps 

M-1 

3.23E-08 

7.43E-02 

3.32E-08 

5.22E-08 

8.74E-08 

3.99E-08 

8.45E-07 

2.41E-06 

3.98E-06 

ps 

MO 

4.33E-09 

1.37E-08 

7.43E-02 

3.31E-08 

1.13E-07 

8.15E-08 

3.44E-07 

5.07E-07 

2.47E-06 

ps 

Ml 

9.92E-09 

2.10E-09 

1.37E-08 

7.43E-02 

1.03E-07 

1.84E-07 

3.66E-07 

9.32E-08 

1.40E-06 

PLi 

7.00E-07 

1.19E-06 

2.18E-06 

3.50E-06 

9.21E-01 

l.OlE-04 

1.07E-04 

5.55E-05 

9.18E-04 

^-1 

2.61E-06 

9.96E-08 

3.85E-07 

2.15E-06 

2.58E-05 

9.21E-01 

5.22E-05 

1.14E-04 

5.55E-05 

ps 

Mo 

8.93E-07 

3.21E-07 

3.05E-07 

7.39E-07 

4.20E-06 

8.07E-06 

9.21E-01 

5.22E-05 

1.07E-04 

ps 

Ml 

1.63E-07 

1.64E-07 

6.13E-08 

4.53E-08 

5.20E-07 

2.14E-06 

8.07E-06 

9.21E-01 

l.OlE-04 

ps 

M2 

3.73E-08 

9.21E-08 

9.97E-08 

7.33E-08 

1.80E-06 

5.20E-07 

4.20E-06 

2.58E-05 

9.21E-01 


Table 4.5: The magnitude of T-matrix elements of a spherical thin-shell transducer computed in 
COMSOL 


matrix can be determined by examining the non-diagonal elements because the applied single¬ 
component of incident wave should theoretically result in only the same harmonic component 
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of the scattered wave. Therefore, the diagonal values of T-matrix should be non-zero for a 
spherical thin-shell transducer, while the off-diagonal values should vanish [12,14]. 

To investigate the off-diagonal elements accurately, the magnitude of the off-diagonal elements 
in the T-matrix were normalized by the magnitude of diagonal elements in the same column 
and the result is shown in Table 4.6. The magnitude of the diagonal elements are significantly 
greater than non-diagonal elements as expected. The largest “leakage” is l.lOE-04 and occurs 


Normalized T-matrix 



pi 

Mo 

pi 

M-1 

pi 

MO 

pi 

Ml 

pi 

M-2 

pi 

M-1 

pi 

Mo 

pi 

Ml 

pi 

M 2 

Pm 

1 

3.50E-07 

1.58E-07 

9.92E-07 

8.15E-08 

4.90E-07 

2.10E-06 

7.59E-06 

3.05E-06 

Pt-i 

4.67E-07 

1 

4.47E-07 

7.03E-07 

9.48E-08 

4.34E-08 

9.18E-07 

2.62E-06 

4.32E-06 

ps 

MO 

6.26E-08 

1.85E-07 

1 

4.46E-07 

1.23E-07 

8.85E-08 

3.74E-07 

5.51E-07 

2.68E-06 

ps 

Ml 

1.43E-07 

2.82E-08 

1.84E-07 

1 

1.12E-07 

2.00E-07 

3.97E-07 

l.OlE-07 

1.52E-06 

P2-2 

l.OlE-05 

1.60E-05 

2.93E-05 

4.72E-05 

1 

l.lOE-04 

1.17E-04 

6.02E-05 

9.96E-04 

ps 

M-1 

3.78E-05 

1.34E-06 

5.18E-06 

2.90E-05 

2.80E-05 

1 

5.66E-05 

1.24E-04 

6.02E-05 

ps 

Mo 

1.29E-05 

4.33E-06 

4.11E-06 

9.95E-06 

4.56E-06 

8.76E-06 

1 

5.66E-05 

1.17E-04 

ps 

Ml 

2.36E-06 

2.21E-06 

8.25E-07 

6.10E-07 

5.64E-07 

2.33E-06 

8.76E-06 

1 

l.lOE-04 

ps 

M 2 

5.40E-07 

1.24E-06 

1.34E-06 

9.87E-07 

1.96E-06 

5.64E-07 

4.56E-06 

2.80E-05 

1 


Table 4.6: The magnitude of the non-diagonal elements in the T-matrix were normalized to the 
magnitude of diagonal elements in the same column 


with quadrupole-quadrupole spherical harmonic component interactions. It is remarkable that 
the leakages are so small and demonstrate the quality of the T-matrix calculations. 

The diagonal elements of the T-matrix are shown in Table 4.7. For the spherical thin-shells. 


ELEMENT 

REAL PART 

IMAG PART 

MAGNITUDE 

PHASE(DEG) 

Tn 

-4.778E-03 

-6.895E-02 

6.911E-02 

-9.40E-I-01 

T22 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E-I-01 

T33 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E-I-01 

T 44 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E-I-01 

T55 

-8.462E-01 

-3.642E-01 

9.212E-01 

-1.57E-I-02 

Tf >6 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E-I-02 

Til 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E-I-02 

Tss 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E-I-02 

T99 

-8.462E-01 

-3.642E-01 

9.212E-01 

-1.57E-I-02 


Table 4.7: Diagonal T-matrix elements for a spherical thin-shell transducer 


the off-diagonal elements of the T-matrix computed with COMSOL were fairly close to zero. 
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Therefore, the off-diagonal elements were neglected in array calculations and the T-matrix for 
the spherical thin-shell transducer is expressed as 


Tsphere “ 


-0.005-0.069i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006 + 0.0741 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006 + 0.0741 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006 + 0.0741 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846- 

-0.364i 0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846-0.364i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846-0.364i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846-0.364i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846 


4.3 The Flextensional Ring-Shell Transducer 

The second major objective of this thesis is to model the performance of a linear array of flex¬ 
tensional ring-shell transducer elements. As no analytical solution is available, an accurate 3D 
model must be constructed to calculate the T-matrix elements of the class V flextensional ring- 
shell transducer in COMSOL. The excellent result obtained for the thin shell FEM calculations 
lent confidence that high accuracy results could be obtained for the ring shell transducer as well. 
Results obtained with COMSOL were compared to prior work as well as to the manufacturer’s 
data in order to investigate the accuracy of the model [7,15]. 

Previous research is provided in the thesis work of Rogerio Pinto [7]. It should be noticed that 
an ultimate goal of this research is to model the performance of the “double-billboard“ array of 
flextensional transducers that was built and tested circa 1990 [1]. 


4.3.1 Geometry 

The dimensional information of the flextensional ring-shell transducer differs slightly from that 
of prior research. The dimensional nodes from the previous attempt to model the 3D flex¬ 
tensional ring-shell in ATILA FEM were used as a reference for the present geometry. The 
coordinate values of the nodes were imported to Microsoft Excel and the required dimensions 
of the geometry were then computed. 

The elements were built using 3D geometric primitives in COMSOL (e.g., sphere, circle, etc.) in 
order to simplify the geometry. The ring-shell transducer was located inside of a fluid spherical 
domain with radiating boundaries at a total radius of 2.5 m. 

The dimensions of the geometry are shown in Table 4.8. 


40 



Parameter 

Value [cm] 

Shell outer radius 

67.32 

Shell height from symmetry center 

14.3 

Shell thickness 

1.04 

Piezoelectric ring edge 

8.0 

Inner radius of piezoelectric ring 

31.9 

Gutter radius of piezoelectric ring 

39.9 

Fiber wrap height 

8.0 

Fiber wrap width 

0.5 

Radius of water domain 

250 


Table 4.8: The ring-shell transducer dimensional parameters 


4.3.2 Material 

Piezoelectric ceramic ring 

The piezoelectric ceramic ring has 144 tangentially-poled lead zirconate titanate ceramic plates, 
separated by 72 steel wedges. The modeling for each tangentially-poled ceramic plate requires 
coordinate system orientation (Figure 3.2) and this process increases the computation time. 
Therefore, “smeared” material properties were used for the composite piezoelectric ceramic 
ring and steel wedges in simplifying the model. Also, the polarization direction altered from 
the tangential direction into the axial direction. 

The properties of the smeared piezo-ring, shown in Figure 4.7 were obtained by McMahon 
and Armstrong through in-air measurements of the segmented ring during various stages of 
assembly [7]. The effect of compression due to fiberglass wrapping and material losses were 
included in the piezoelectric material properties. 



Figure 4.7: Smeared piezo-ring 
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A further transformation is required of the elasto-piezoelectric material property matrix in order 
to model a circumferentially-polarized ring as an axially-polarized ring. The required steps are 
explained in Appendix A. After transformation, the elasto-piezoelectric matrix becomes 


Si ■ 


r 

■^11 

■^13 

il2 

0 

0 

0 

0 

0 

d3i 


■ Ti ■ 

S2 


■^13 

^33 

■^13 

0 

0 

0 

0 

0 

d33 


Ti 

S3 


5i2 

■^13 

■^ll 

0 

0 

0 

0 

0 

d3\ 


T3 

S4 


0 

0 

0 

‘^44 

0 

0 

^15 

0 

0 


T4 

Ss 

= 

0 

0 

0 

0 

‘'66 

0 

0 

0 

0 


Ts 

S6 


0 

0 

0 

0 

0 

S 44 

0 

di5 

0 


Te 

Di 


0 

0 

0 

^15 

0 

0 

*=^11 

0 

0 


El 

Di 


0 

0 

0 

0 

0 

^15 

0 

^11 

0 


El 

D 3 _ 


_ <^31 

d33 

d3i 

0 

0 

0 

0 

0 

^33 J 


.E3 . 


(4.8) 


for which the smeared material properties are: 


Parameter Value Unit Notation 


Density 

Compliance Matrix 

7550 

kg/m^ 

.E 

^11 

8.202-0.32808i 

lO-'^mVN 

^12 

-2.703-1-0.10812i 



-3.538-1-0.141521 


.E 

^33 

11.99-0.47961 

lO-'^mVN 

.E 

^44 

26-1.04i 

lO-'^mVN 

.E 

^66 

21.8-8.72i 

lO-'^mVN 

Coupling Matrix 

d3i 

-1 

10-‘®C/N 

d33 

2.139 

10-'®C/N 

d\5 

-2.811 

10-'®C/N 

Relative Permittivity 

efi/eo 

1297-5.1861 


4/^ 

1297-5.1861 


£0 

8.85 

lO-'^F/M 


Table 4.9: The flextensional transducer smeared piezoelectric material properties 


To select an axial polarization direction the user-defined base vector system was added in the 


42 



definition branch of the model and the cylindrical coordinate system was defined using the unit 
basis vectors (i, j, k) [21] 

xi = = cos + sin 

*2 = = — sin + cos 0j (4.9) 

X3 = Bz = k . 

It should be noted that this modification of polarization requires adjustment of the computed 
pressure values. The required correction was described in the masters thesis of Rogerio Pinto 
as the following step [7] 

“divide the displacement and pressure fields by the ratio of the circumferential length of a 
‘smeared’ piezoelectric element (which is equal to 1/144 of the the ring circumferential length) 
to its height. In the present case this ratio is 0.1958.” 

The correction factor calculated using this ratio as -14.16 dB, i.e. a value of 14.6 dB was added 
to the computed pressures, for a given voltage drive. Note this correction is only required for 
computing the radiated pressure field; it does not apply for the short-circuit scattering calcula¬ 
tions required to compute the T-matrix values. 

Shell 

The shell, shown in Figure 4.8 was made of a special alloy steel (ST4340) and modeled by the 
user defined material properties given in Table 4.10, and taken from an input data file of ATILA 
FEM [7]. 



Figure 4.8: ST4340 shell 


Fiberglass 

The fiberglass sheath, shown in Figure 4.9, provides compression to the ring, having been ap¬ 
plied under tension and bonded with epoxy resin, and the smeared piezoelectric ceramic proper- 
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Parameter 

Value 

Notation 

Density 

7700 kg/m^ 

P 

Young’s Modulus 

195 GPa 

Y 

Poisson’s ratio 

0.28 

V 


Table 4.10: Shell material properties 


ties include the effect of the fiberglass wrapping compression. For that reason, fiberglass wrap¬ 
ping is added without any specified compression to the model; Only its linear elasto-dynamic 
properties are specified, as shown in Table 4.11. 


Parameter 

Value 

Notation 

Density 

2000 kg/m^ 

P 

Young’s Modulus 

10 GPa 

Y 

Poisson’s ratio 

0.28 

V 


Table 4.11: Fiberglass material properties 
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Fluid Domain 

The fluid domain has the same material properties as given in Table 4.2. The fluid media include 
exterior water and interior air. 

4.3.3 Meshing 

The maximum element size is defined with respect to the Nyquist criteria which states that the 
maximum element size must be smaller than A/5 (This criteria was applied with respect to 
the speed of sound in water). The frequency used is 605 Hz which is close to the operating 
frequency and also provides “^ 0 = 1 ”. 

The upper surface of the piezo-ring was meshed with free triangular mesh elements using the 
size adjustment of “extra fine.” The meshed surface was swept by a built-in feature of COMSOL 
through the axial direction to the bottom surface of the piezo-ring. 



(a) (b) 

Figure 4.10: (a) Converted faces with quadrilateral mesh (b) Swept mesh with free triangular 
elements 

The number of the elements swept in were adjusted to be 10 equal hexahedron elements using 
the “distribution” node feature of COMSOL (Figure 4.10b). The surface boundary between ring 
faces and the surrounding fluid was converted to a quadrilateral mesh to generate an adjacent- 
free tetrahedral mesh for the surrounding domain (Figure 4.10a). 



Figure 4.11: Complete mesh of the ring-shell transducer 
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The remaining geometry was meshed using free triangular mesh elements with a size adjustment 
of “finer” (Figure 4.11). The complete mesh, shown in Figure 4.12 consists of 207,003 elements 
and the average element quality of the entire geometry is 0.66. 



Figure 4.12: Complete mesh of entire geometry 


Geometry Entity Level 

Element Type 

Amount 

Avg. Element Quality 

All Domain 

Tetrahedral 

200,483 

0.67 


Pyramid 

2,680 

0.20 


Prism 

3,840 

0.48 

Boundaries 




Radiation 

Triangular 

5,056 

0.99 

Piezo ceramic ring 

Triangular 

3,496 

0.54 

Wet surface of transducer 

Triangular 

7,588 

0.91 

Entire Geometry 


207,003 

0.66 


Table 4.12: Meshing properties of the 3D spherical transducer model with extra fine mesh for 
the piezo-ring and with a size adjustment of finer for the remaining geometry 


4.3.4 Validation of Model 

To quantify validation of the model, the source level (SL) at 3000 Vrms at 1 m was computed 
in COMSOL, and the correction factor, -14.16 dB, was added to the computed SL (See Chapter 
4.3.2 for details). 
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Figure 4.13 shows the SL as a function of the frequency computed by COMSOL using the 
baseline properties at 10 Hz intervals, from 200 Hz to 2000 Hz and Figure 4.14 shows a com¬ 
parison of the SL computed in ATILA FEM [7] against the manufacturer’s data of the SL as a 
function of frequency. The number of DOFs solved was 495,594 and the solution time for each 
frequency was 255 s. 

The first peak value of 211.6 dB re 1 /tPa/V at 1 m occurred at 605 Hz with 3000 Vrms driven 
voltage. While this value differs by 1.4 dB in SL from the manufacturer’s data (213 dB), it is 
much better than ATILA FEM results which differed by 4 dB from the measured value at the 
resonance frequency. 


SOURCE LEVEL (SL) DRIVEN AT SOOOVrms 



Figure 4.13: SL at Im on acoustic axis driven 3000 Vrms, obtained with COMSOL ACPZ 
interface (main model). 


The very narrowband irregularities observed in Figure 4.13 (660 Hz, ~1200 Hz, ~ 1700 Hz ) 
were found to be uncontrolled shell modes. They are not seen in experimental measurements. 

Figure 4.15 depicts the bandwidth at resonance frequency, and it shows the interpolated value 
of ±3dB with respect to computed values. The spurious value at 660 Hz was omitted as being 
misleading. The mechanical quality factor, Qm, was calculated 4.35 and a bandwidth of 139 
Hz. The Qm differ 0.4 and the bandwidth 16 Hz from the manufacturer’s measured value. 

The model shows (Figure 4.13) the second resonance of the developed model occurring at 950 
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Figure 4.14: SL on acoustic axis driven 3000 Vrms, obtained from the master thesis of 
Pinto, 1993 [7] 


TRANSMITTING VOLTAGE RESPONSE(TVR) DRIVEN AT 3000Vrms 



Figure 4.15: The bandwidth at primary resonance frequency(main model) 


Hz which is in excellent agreement with manufacturer’s data, but the model displays an approx¬ 
imately 6 dB higher peak after the secondary resonance frequency. As the frequency increases 
to 1000 Hz, the model shows some discrepancies. After 1600 Hz the model shows a lack of 
convergence to manufacturer’s data. 
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At first glance, the expected reason for this deviation may be due to internal material losses. The 
internal material losses were not included in the model and could vary the computed SL values. 
Secondly, the mesh quality (Figure 4.16) may affect the results at higher frequencies since 
meshing quality is related to operating wavelength. The proper meshing size at low frequency 
becomes too coarse at high frequency. Therefore, increasing the frequency, reduces the mesh 
quality of the model. It should be noted, that the mesh size of the developed model was adjusted 
according to “/=474 Hz” and the selected mesh size was not refined enough for high frequency. 

Further investigation of meshing quality indicated another issue: the meshing of thin structures 
deteriorated. Coupling of a thin structure with a larger object obstructs the growth rate of the 
mesh size. A sample of this situation occurs with the coupling between fiberglass wrapping and 
the piezo-ring domain. The quality distribution of the mesh is shown in Figure 4.16. 


0.2 

0.1 

0 

▼ 0.0155 

Figure 4.16: Mesh quality plot of the main model for the ring-shell transducer. 

Although, the model shows some deviation above the secondary resonance frequency, the model 
still demonstrates quite good agreement with the manufacturer’s data between 400-1000 Hz 
(“A:a=l”) and better agreement than the previous model (Pinto,1993). 

The developed model is sufficient for the purpose of a T-matrix calculation because the perfor¬ 
mance of the class V flextensional transducer will reside in vicinity of “A:a=l”, and so was used 
as the base model for the T-matrix calculation. 

4.3.5 T-Matrix Calculation 

The spherical harmonic components of an incident standing wave were applied at the fluid/solid 
surface of the transducer, and the spherical harmonic scattered amplitudes were determined 
from the resultant scattered pressure field in the fluid. 

Post processing properties were used to create a spherical parametrized surface within the water 
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domain at radii of 0.5 m. Figure 4.17 shows the defined spherical parametrized surface in water 
domain. 


The magnitude of the T-matrix elements are provided in Table 4.13. 



Figure 4.17: A parametrized spherical surface was created at radii of 0.5 m (shown in red; The 
great circles are on the radiation boundary at 2.5 m). 


T-matrix Elements 



pi 

MO 

pi 

M-l 

pi 

MO 

pi 

Ml 

pi 

M-2 

pi 

M-l 

pi 

MO 

pi 

Ml 

pi 

M2 

ps 

^00 

9.30E-01 

2.76E-05 

1.64E-04 

7.68E-05 

4.93E-08 

2.63E-06 

1.93E-02 

1.64E-05 

2.25E-05 

ps 

M-1 

7.54E-05 

5.06E-02 

1.63E-05 

6.84E-05 

4.62E-07 

3.03E-06 

1.17E-06 

3.73E-06 

3.93E-06 

ps 

MO 

5.07E-04 

1.74E-05 

5.68E-02 

6.19E-05 

2.79E-07 

1.38E-06 

1.23E-05 

1.24E-05 

7.02E-06 

ps 

Ml 

8.49E-05 

1.34E-06 

1.71E-05 

7.55E-02 

9.25E-08 

1.87E-07 

2.16E-06 

6.73E-06 

1.86E-06 

ps 

M-2 

1.33E-04 

2.67E-06 

l.llE-05 

5.63E-05 

4.90E-03 

1.75E-07 

3.37E-06 

5.95E-06 

6.81E-05 

ps 

M-1 

1.58E-05 

1.22E-05 

2.21E-05 

7.87E-06 

2.21E-07 

4.45E-03 

1.77E-06 

1.56E-06 

6.53E-06 

ps 

MO 

9.57E-02 

3.36E-06 

2.03E-05 

9.41E-06 

7.42E-09 

3.12E-07 

4.23E-03 

1.74E-06 

2.70E-06 

ps 

Ml 

6.35E-06 

3.42E-07 

2.29E-06 

8.47E-06 

3.68E-08 

8.86E-08 

1.21E-07 

4.45E-03 

8.64E-08 

ps 

M2 

2.00E-06 

1.32E-07 

3.86E-07 

1.46E-06 

5.84E-08 

3.78E-08 

4.40E-08 

2.00E-07 

4.90E-03 


Table 4.13: The magnitude of T-matrix elements of a class V flextensional ring-shell transducer 
computed in COMSOL 


For convenience of analysis, the magnitude of the non-diagonal elements in the T-matrix were 
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normalized to the magnitude of diagonal elements in the same column and the results are shown 
in Table 4.14. 


Normalized T-matrix 



pi 

MO 

pi 

M-1 

pi 

MO 

pi 

Ml 

pi 

M-2 

pi 

M-1 

pi 

MO 

pi 

Ml 

pi 

M2 

ps 

^00 

1 

5.46E-04 

2.88E-03 

1.02E-03 

l.OlE-05 

5.91E-04 

1 

3.68E-03 

4.59E-03 

ps 

M-1 

8.11E-05 

1 

2.87E-04 

9.06E-04 

9.42E-05 

6.81E-04 

6.06E-05 

8.39E-04 

8.01E-04 

ps 

MO 

5.45E-04 

3.44E-04 

1 

8.20E-04 

5.68E-05 

3.09E-04 

6.38E-04 

2.78E-03 

1.43E-03 

ps 

Ml 

9.13E-05 

2.65E-05 

3.02E-04 

1 

1.89E-05 

4.20E-05 

1.12E-04 

1.51E-03 

3.80E-04 

ps 

M-2 

1.43E-04 

5.28E-05 

1.96E-04 

7.46E-04 

1 

3.93E-05 

1.74E-04 

1.34E-03 

1.39E-02 

ps 

M-1 

1.70E-05 

2.42E-04 

3.89E-04 

1.04E-04 

4.51E-05 

1 

9.17E-05 

3.50E-04 

1.33E-03 

ps 

MO 

1.03E-01 

6.64E-05 

3.58E-04 

1.25E-04 

1.51E-06 

7.02E-05 

2.19E-01 

3.91E-04 

5.50E-04 

ps 

Ml 

6.83E-06 

6.76E-06 

4.03E-05 

1.12E-04 

7.50E-06 

1.99E-05 

6.29E-06 

1 

1.76E-05 

ps 

M2 

2.16E-06 

2.60E-06 

6.79E-06 

1.94E-05 

1.19E-05 

8.50E-06 

2.28E-06 

4.48E-05 

1 


Table 4.14: The magnitude of non-diagonal elements in the T-matrix normalized to the magni¬ 
tude of diagonal element in the same column 


In the vicinity of the primary resonance frequency, 605 Hz, the flextensional transducer radiates 
both quadrupole and monopole pressure fields. It is found that the linear quadrupole radiated 
pressure coefficient is about 1/10 of the monopole coefficient ( see section 5.2 for details) When 
a monopole field is incident, the ring-shell generates a quadrupole field and vice versa, charac¬ 
teristic of the resonant vibration. This general dynamic response of the transducer can be seen 
from the T-matrix table. It is remarkable that the magnitude of Tjj (quadrupole-quadrupole) is 
smaller than the magnitude of T\q (quadrupole-monopole) which are highlighted in Table 4.14. 
Thus, the calculated T-matrix captures the interaction between monopole and quadrupole exci¬ 
tation. 

The T-matrix for the ring-shell transducer is not diagonal while the T-matrix of the spherical 
thin-shell transducer essentially is. This is expected behavior as the geometry of the flexten¬ 
sional transducer is not spherical. It should be noted, that the T-matrix is a kind of transfer 
function. The diagonal elements are represented in Table 4.15. 

Similar to the spherical thin-shell transducer, the T-matrix of the flextensional ring-shell trans¬ 
ducer is nearly a diagonal matrix except for elements of monopole-linear quadrupole interac¬ 
tions (ri 7 ,r 7 i). Unlike the spherical thin-shell transducers, the monopole-linear quadrupole 
term must be included in the T-matrix. Neglecting small values in the matrix, the T-matrix of 
the flextensional ring-shell becomes. 
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ELEMENT 

REAL PART 

IMAG PART 

MAGNITUDE 

PHASE(DEG) 


7^11 

-9.237E-01 

-1.088E-01 

9.301E-01 

-1.73E-I-02 

Til 

-1.538E-02 

-4.820E-02 

5.060E-02 

-1.08E-I-02 

T 33 

-1.652E-02 

-5.435E-02 

5.680E-02 

-1.07E-I-02 

T 44 

-2.467E-02 

-7.134E-02 

7.548E-02 

-1.09E-I-02 

T 55 

-1.674E-03 

-4.607E-03 

4.902E-03 

-l.lOE-l-02 

766 

-1.660E-03 

-4.129E-03 

4.450E-03 

-1.12E-I-02 

Til 

-3.586E-03 

-2.242E-03 

4.229E-03 

-1.48E-I-02 

Tss 

-1.660E-03 

-4.129E-03 

4.450E-03 

-1.12E-I-02 

T 99 

-1.674E-03 

-4.607E-03 

4.902E-03 

-l.lOE-i-02 


Table 4.15: Diagonal T-matrix elements for a class V flextensional ring-shell transducer 


TRing shell — 


-0.924-0.109i 0 0 0 0 0 -0.019 - 0.002i 0 0 

0 -0.015-0.048i 0 0 0 0 0 0 0 

0 0 -0.017-0.054i 0 0 0 0 0 0 

0 0 0 -0.025-O.OVli 0 0 0 0 0 

0 0 0 0 -0.002-0.005i 0 0 0 0 

0 0 0 0 0 -0.002 - 0.004i 0 0 0 

-0.095 - 0.009i 0 0 0 0 0 -0.004 - 0.002i 0 0 

0 0 0 0 0 0 0 -0.002-0.004i 0 

00000000 -0.002 - 0.005i 


4.4 Improvement on the Ring-Shell Transducer Model at Higher 
Operating Frequencies 

Although the accuracy of the developed model was verified up to 1000 Hz, corresponding to 
the frequency of interest in this thesis to calculate the T-matrix elements, the SL of the devel¬ 
oped model demonstrated some deviation from experimental values at higher frequencies (Fig¬ 
ure 4.13). To achieve more accurate results at higher frequencies, improvements were made 
to the first developed model. The first model developed is referred to as the “main model” 
throughout this section. 

4.4.1 Refined Mesh Model 

For frequencies up to 2 kHz, the wavelength of sound in water is 0.745 m (c=1490 m/s). The 
maximum meshing element size should be specified with respect to X/5 at this frequency. 

In this limit, the maximum element size of meshed geometry equals 0.149 m which is much 
smaller than the maximum mesh element size of the main model (0.275 m). To obtain an ac¬ 
curate result at higher frequency, the mesh of the geometry was refined and adjusted satisfy the 
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A/5 criteria at 2000 Hz. The refined mesh consist of 413,778 elements and 733,161 DOFs. 

The SL of the refined mesh model is shown in Figure 4.18 and the SL was ealeulated at 10 Hz 
intervals, from 200 Hz to 2000 Hz. Primary and secondary resonances were captured at the 
same frequency as those in the main model, as are the amplitudes. Mesh refinement did not 
produee any signifieant improvement. Thus, the general behavior of SL stayed the same and the 
refined mesh eould not remove diserepancies between the graph of computed SL values and the 
graph of measured SL values (Figure 4.14). It should also be noted that the computation time 
increased due to the inereased number of DOFs of the refined mesh. 



Figure 4.18: SL at 1 m on acoustic axis is driven at 3000 Vrms of the refined mesh model, the 
mesh size is adjusted with respect to 2000 Hz. 

4.4.2 Improved Model 

Upon further investigation, it was revealed that the mesh quality was reduced in the neigh¬ 
borhood of the fiberglass wrapping and the surface of the shell. The fiberglass wrapping was 
re-meshed with triangular elements and swept through the axially direction from top to bottom 
with a limited maximum mesh element size to A/30 or 0.025 m at 2500 Hz. 

In addition to mesh improvement, the spherical wave radiation boundary condition was replaced 
with a spherieal perfeetly matehed layer (PML) to absorb the ineident waves without refleetion 
as they strike the radiation boundary. The PML surrounded the fluid domain and had a 0.5 m 
thickness, reducing the radius of fluid domain to 2 m. The updated geometry is represented in 
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Figure 4.19. The PML mesh was swept radially with a distribution of ten hexahedron elements 
through its thickness, and the water domain between PML and ring-shell was meshed with a 
boundary layer mesh. The boundary layer mesh is recommended by COMSOL in order to 
enhance the precision of the far field calculation near the PML [17]. The complete mesh consist 
of 238,503 elements and 792,888 DOFs. 



PML 


NG-SHELL 


Figure 4.19: The geometry of the improved model, PML surrounded the fluid domain. 


Figure 4.20 shows the SL of the improved model where again the SL was calculated at 10 
Hz intervals, from 200 Hz to 2000 Hz. The behavior of the model changed. The secondary 
resonance frequency shifted left by 35 Hz from the manufacturer’s data to 985 Hz, and the SL 
became higher than the manufacturer’s data by 6 dB. 

The effects of the improvements were seen at frequencies higher than the secondary resonance. 
Beyond that frequency, the SL agrees with the manufacturer’s data (Figure 4.14) and the dif¬ 
ference was reduced to 2 dB, showing fairly good agreement with the manufacturer’s data up 
to 1600 Hz. For frequencies higher than 1600 Hz, the SL again deviates slightly from the 
manufacturer’s data. 
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SOURCE LEVEL (SL) DRIVEN AT 3000Vrms 



Figure 4.20: SL on acoustic axis driven at 3000 Vrms of the improved model, PML and bound¬ 
ary layer mesh added. 


4.4.3 Result 

In addition to the main model, two different models were examined for the purpose of improving 
accuracy at higher frequencies. Figure 4.21 depicts the comparison of SL of the main model, 
the refined mesh model and the improved model. 

The main model and refined mesh models performed identically from 400 Hz to 800 Hz. The 
improved model produced a lower SL than the other models at higher frequencies, and it con¬ 
verged with the manufacturer’s data. Thus, the improved model performed better when com¬ 
pared to the other model, although it captured the secondary resonance with a 35 Hz difference 
from experiment. 

Reasons for the improvement in the SL value are due in part to the changes in the mesh structure 
and applied PML as boundary condition. 


SOURCE LEVEL (SL) DRIVEN AT 3000Vrms 



Figure 4.21: The SL comparison of three developed model 
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The detailed SL of frequency interest in the neighborhood of the first resonance is illustrated 
in Figures 4.23 through 4.24. The comparison of the bandwidth and the Q values are shown in 
Table 4.16. The three models performed nearly the same and show fairly good agreement with 
the manufacturer’s data between 250-1000 Hz. 


Model 

BW (Hz) 

BW Error(Hz) 

BW(%) 

Q 

Q Error 

Q Error(%) 

Main 

138.5 

16.5 

10.64 

4.33 

0.43 

11 

Refined 

139.5 

15.5 

10.00 

4.30 

0.4 

10.26 

Improved 

139 

16 

10.03 

4.35 

0.45 

10.35 


Table 4.16: comparison of bandwidth and mechanical Q value from three model against manu¬ 
facturer’s data 


All the models reproduced equally well the SL up to 1000 Hz. For this reason, the main model 
was preferred for the T-matrix element calculation because of reduced computation time for 
equivalent accuracy. 


Model 

number of element 

average 

worst 

DOFs 

Computation Time 

main 

207,003 

0.66 

2.29E-05 

495594 

258 s. 

refined 

413,778 

0.73 

4.23E-05 

733161 

353 s. 

improved 

238,503 

0.56 

2.23E-06 

792888 

550 s. 


Table 4.17: The computation time values and the general mesh information of the three model 


SL DRIVEN AT 3000Vrms(MAIN MODEL) 



Figure 4.22: The bandwidth at primary resonance frequency (main model) 
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Figure 4.23: The bandwidth at primary resonance frequency (refined model) 


SL DRIVEN AT SOOOVrms (IMPROVED MODEL) 



Figure 4.24: The bandwidth at primary resonance frequency (improved model) 
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CHAPTER 5: 

Linear Array Performance Modeling by Using T-Matrix 

Method 


The vertical linear array with the spherical thin-shell transducer elements and the vertical linear 
array with the ring-shell transducer elements were programmed in MATLAB, using the T-matrix 
method. The required T-matrix elements for both types of the transducers were reported in 
Chapter 4. 

5.1 The Development Sequence of an Array Model 

1. Computation of the free-held radiated pressure coefficients 

2. Calculation of the translation matrix using the addition theorem 

3. Implementation of the T-matrix to calculate the scattered pressure coefficients 

4. Computation of the total acoustic pressure held and generation of the far-held beam pat¬ 
tern 


5.2 Computation of the Free-Field Radiated Pressure Coefficients 


The free-held radiated pressure coefficients of both transducer types were calculated at a 1 V 
peak driven potential. That a driving peak was used versus 1 Vrms is of no consequence as 
only normalized beam patterns are presented in this chapter. The free-held radiated pressure 
amplitudes were computed by solving equation (2.40) using COMSOL. The number of spher¬ 
ical harmonic components kept was truncated at the the quadrupolar degree, n-2, so the total 
number of spherical harmonic components kept was 9 (n=0,l,2 and m--n..n). 

Results for the spherical thin-shell transducer are listed in Table 5.1, and the results for the 
ring-shell transducer are given in Table 5.2. 
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operation freq 

474Hz(ka=l) 





Mode 

Real 

Imaginary 

Amplitude 

Amplitude 

IPool 

Angle(rad) 

pr 

^00 

1.349E+00 

1.946E-I-01 

1.951E-I-01 

l.OOOE+00 

1.502E-I-00 

pr 

M-1 

6.323E-06 

6.863E-06 

9.331E-06 

4.784E-07 

8.263E-01 

pr 

MO 

3.608E-07 

1.330E-06 

1.378E-06 

7.067E-08 

1.306E-h00 

pr 

Ml 

2.909E-06 

-9.429E-07 

3.058E-06 

1.568E-07 

-3.135E-01 

pr 

M-2 

2.633E-04 

1.202E-04 

2.894E-04 

1.484E-05 

4.285E-01 

pr 

M-1 

-6.943E-04 

3.410E-04 

7.735E-04 

3.966E-05 

2.685E-h00 

pr 

MO 

2.027E-04 

9.465E-05 

2.237E-04 

1.147E-05 

4.369E-01 

pr 

Ml 

3.074E-05 

3.885E-05 

4.954E-05 

2.540E-06 

9.015E-01 

pr 

M2 

-1.306E-06 

-8.569E-06 

8.668E-06 

4.444E-07 

-1.722E+00 


Table 5.1: The radiated pressure coefficients of the spherical thin-shell transducer for monopole, 
dipole, and quadrupole {n-0, 1, 2 and m=-n..n), for 1 V peak drive. 


operation freq 

605Hz(ka=l) 





Mode 

Real 

Imaginary 

Amplitude 

Amplitude 

|P|iol 

Angle (rad) 

pr 

MO 

6.431E-I-00 

9.823E-01 

6.506E+00 

l.OOOE+00 

1.516E-01 

pr 

M-1 

1.026E-04 

-6.191E-04 

6.275E-04 

9.645E-05 

-1.407E-I-00 

pr 

MO 

1.109E-03 

3.419E-03 

3.594E-03 

5.524E-04 

1.257E-I-00 

pr 

Ml 

4.475E-04 

-4.552E-04 

6.384E-04 

9.812E-05 

-7.939E-01 

pr 

M-2 

-4.804E-06 

1.083E-03 

1.083E-03 

1.665E-04 

1.575E-I-00 

pr 

M-1 

1.188E-04 

1.791E-04 

2.149E-04 

3.303E-05 

9.853E-01 

pr 

MO 

6.738E-01 

3.858E-02 

6.749E-01 

1.037E-01 

5.719E-02 

pr 

Ml 

-2.213E-05 

2.659E-05 

3.459E-05 

5.317E-06 

2.265E-I-00 

pr 

M2 

8.583E-06 

-6.395E-06 

1.070E-05 

1.645E-06 

-6.403E-01 


Table 5.2: The radiated pressure coefficients of the ring-shell transducer for monopole, dipole, 
and quadrupole {n-0, 1, 2 and m--n..n), for 1 V peak drive. 


Table 5.1 indicates that the free-field radiated pressure amplitudes of the spherical thin-shell 
transducer are negligible except for the monopole expected. Because of symmetry, 

coefficients other than Pqq are close to zero, indicating the accuracy of the finite element model. 

Similarly, Table 5.2 shows that the free-field radiated pressure amplitudes of the ring-shell trans¬ 
ducer are negligible except for the monopole (Pqq) and linear quadrupole (PIq), the monopole 
coefficient being about “10” times the linear quadrupole coefficient. This is expected due to the 
nature of the flexural motion of the ring-shell transducer. 
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Negligible values of the free-field radiation are ignored in far-field calculation. The free-field 
radiated pressure vectors then become for a 1 V peak drive, 


pR — 

^ Sphere 


1.345+19.459i 
0 
0 
0 
0 
0 
0 
0 
0 


pR 

Ring—shell 


6.431+0.982/ 

0 

0 

0 

0 

0 

0.674 + 0.039/ 
0 
0 


and they were imported to MATLAB for further calculation. 


(5.1) 


(5.2) 


5.3 Calculation of the Translation Matrix 

The translation matrix, “Giq”, is needed to relate the outgoing wave component from 

transducer to the “n, free-field incident component on transducer The translation 
matrix is obtained from a spherical addition theorem (See Chapter 2.4.5 for details). Bold 
capital “G” denotes a translation matrix from one transducer to another in the array 
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G/,00 G/,00 G;,oo G/,00 G;,oo G/,oo G;,oo G/,oo G/,oo 

^,00 ^,1 — 1 ^,10 ^,11 2 ^, 2—1 ^,20 ^,21 ^,22 

G/,1-1 G/,1-1 G/,1-1 Gz,i-i G/,1-1 G/,1-1 G/,1-1 Gi,i-i Gi,i-i 

q,00 q,l-l q,l0 q,n q,2-2 q,2-l q,20 q,2l q,22 

Gz,io G/,10 G;,io G/,10 G;,io G/,io G;,10 G/,io G;,10 

^,00 q,l-l q,l0 q,ll q,2-2 q,2-\ q,20 q,2\ q,22 

Gi,n G/,11 G;,11 G/,11 G;,11 G/,ii Gz,ii G;,11 Gz,ii 

^,00 ^,1 — 1 ^,10 q,ll ^)2—2 ^,2—1 ^,20 q,2l q,22 

Gi,2-2 G/,2-2 G/,2-2 G/,2-2 G/,2-2 G/,2-2 G/,2-2 Gi,2-2 G/,2-2 

^,00 ^,1 — 1 ^,10 ^,11 ^)2—2 ^,2—1 ^,20 ^,21 ^,22 

Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 Gl,2-1 

q,00 ^,1 — 1 ^,10 ^,11 ^>2—2 ^,2—1 ^,20 ^,21 ^,22 

G/,20 G/,20 G/,20 G/,20 G/,20 G/,20 G/,20 G/,20 G/,20 

^,00 ^,1-1 ^,10 q,ll q,2-2 q,2-l q,20 q,21 q,22 

Gl,21 G/,21 G;,21 G/,21 G/,21 G/,21 G;,21 G/,21 G/,21 

^,00 q,l-l q,10 q,ll q,2-2 q,2-l q,20 q,21 q,22 

G/,22 G/,22 G/,22 G/,22 G/,22 G/,22 G/,22 G/,22 G/,22 

^,00 ^,1 — 1 ^,10 q,ll ^>2—2 ^,2—1 ^,20 q,21 q,22 


Equation (2.43) was employed in MATLAB to construct the translation matrix, Giq, I ^ q . 


The row refers to the interaction of the ‘7^*” transducer with the other transducers in the array. 
Computation time for producing the translation matrix depends on the number of elements 
within an array 

0 Gi2 • • • Gin 
G21 0 •• • G2N 


GnI Gn2 0 


where 


pout = pR + pS . 






5.4 Implementation of the T-matrix 

The scattered pressure amplitudes can be obtained solving (2.49). The product of both sides 
with the inverse matrix of (I — a) gives the scattered pressure vectors. 



where a is defined for an array which is composed of the identical elements, 


(5.4) 
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a = 


/ 0 TGi 2 • • • TGin \ 

TG21 0 TG2N 

\ TGni TGn2 ••• 0 / 

Considering a three-element transducer array configuration, for example, the free-field scattered 
pressure amplitudes become 


pf 

\ 


I 

-TG12 

-TGi 3 ■ 

-1 

0 

TG12 

TGi 3 


/ pf 

\ 

pf 


— 

-TG21 

I 

-TG23 


TG21 

0 

TG23 


pf 


VPf 

) 


-TG31 

-TG32 
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. TG31 

TG32 

0 


VP3 
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In Chapter 4.2.6 the complex coefficients of the T-matrix for the spherical thin-shell transducer 
were calculated as 


Tsphere ~ 


-0.005-0.069i 0 0 0 0 

0 -0.006+ 0.074i 0 0 0 

0 0 -0.006+ 0.074i 0 0 

0 0 0 -0.006 + 0.074i 0 

0 0 0 0 -0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 


0 0 0 0 

0 0 0 0 

0 0 0 0 

0 0 0 0 

.846-0.3641 0 0 0 0 

-0.846-0.364i 0 0 0 

0 -0.846-0.364i 0 0 

0 0 -0.846-0.364i 0 

0 0 0 -0.846-0.364i 


In Chapter 4.3.5 the complex coefficients of T-matrix for the ring-shell transducer were calcu¬ 
lated as 


TRing shell — 


-0.924-0.109i 0 0 0 0 0 -0.019 - 0.002i 0 0 

0 -0.015-0.048i 0 0 0 0 0 0 0 

0 0 -0.017-0.054i 0 0 0 0 0 0 

0 0 0 -0.025-0.071i 0 0 0 0 0 

0 0 0 0 -0.002-0.005i 0 0 0 0 

0 0 0 0 0 -0.002 - 0.004i 0 0 0 

-0.095 - 0.009i 0 0 0 0 0 -0.004 - 0.002i 0 0 

0 0 0 0 0 0 0 -0.002-0.004i 0 

00000000 -0.002 - 0.005i 
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5.5 Computation of the Total Acoustic Pressure Field and 
Far-Field Beam Pattern 

In the array models, a line of Nx transducer sources with the adjacent elements equally spaced 
were analyzed and are as shown in Figure 5.1. The transducers are aligned with the polar axis 
(z axis). The beam patterns of the developed models, therefore, have rotational symmetry with 
respect to the array (polar) axis. Beam patterns are represented by polar coordinates (0,0), 
and the far-field approximation was used with phase difference between elements, due to their 
position. This was taken into account for the calculations of far-field acoustic field and beam 
pattern. 



(a) (b) 

Figure 5.1: Orientation of the transducers; spacing A/2 and "ka = 1" (a) Spherical thin-shell 
transducer (b) Class V ring-shell transducer 


The total acoustic pressure in the far-field was obtained from superposition of the radiated and 
the scattered complex acoustic pressure amplitudes originating from the transducers in the array 
(See Chapter 2.4.7 for details). Equations (5.5) and (5.7) were employed in MATLAB to obtain 
far-field pressure amplitudes.The radiated pressure amplitudes are reported in Section 5.2, and 
the scattered pressure amplitudes can be obtained from equation (5.4). All equations were 
combined in MATLAB, and the code used is presented in Appendix B. 

The far-field pressure is 
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„-ikr Nx ^ Wmax n 

^ Jk-roi ^ ,•(«+!) 


= E PL+PLn (5.5) 


;=1 n=0 m=—n 


where 


k-foi = k{ xqi sin 6 cos (j) +yoi sin0sin0 + zoi cos0) . 


(5.6) 


and Xqi, you zoi ^rc the rectangular coordinates of the origin of transducer 
The normalized far-field beam pattern is 


b = 201ogio 


( |p//(r,0,(^)| \ 


(5.7) 


To understand the array interaction effect on the far-field beam pattern different linear array 
configurations were examined. For all configurations, 3D normalized far-field beam patterns 
were plotted using the T-matrix method. As the horizontal (i.e., perpendicular array axis) beam 
patterns of all configurations were omni-directional, the detailed investigation of the effect of 
the T-matrix method on beam pattern was focused on the vertical beam (plane including array 
axis direction) patterns. The resultant normalized vertical beam patterns were compared against 
those computed for non-interacting simple sources. 


5.6 Figures of the Far Field Beam Patterns 

The 2D normalized vertical far-field beam patterns were plotted in dB as a function of the an¬ 
gle (in degrees). The beam patterns obtained with the T-matrix method are compared to point 
sources on the same plot. For the polar plot, magnitudes lower than “-30 dB” were omit¬ 
ted since values lower than -30 dB are negligible compared to the maximum intensity. The 
“0 dB” reference value was shifted to 30 dB to improve visualization. Similarly, for improving 
the visualization of cartesian plots of the beam patterns, magnitudes which were lower than 
“-50 dB” were omitted and the “0 dB” reference shifted to 50 dB. The angles “90” and “270” 
correspond to the broadside direction while “0” and “180” correspond to the “array axis (end 
fire)” direction. 

Several configurations were examined, to investigate the performance of the array and to deter¬ 
mine array interaction. 
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5.6.1 The Model of the Linear Array With Spherical Thin-Shell Elements 

Three-Element Array Configuration 

Figure 5.2 depicts the 3D normalized far-field beam pattern of a linear array with three evenly- 
spaced, identical, spherical thin-shell transducer elements at “^0=1” (/=474 Hz). The 
inter-element spacing is equal to A/2 where X is the acoustic wavelength. In terms of am¬ 
plitudes, the broadside is the major lobe, and a nodal surface does not exist. Instead of the 
expected conical nodal surface, some degradation is observed between broadside and end fire. 

Figure 5.3 is the normalized vertical beam pattern of the three evenly-spaced, identical, spher¬ 
ical thin-shell transducers. The inter-element spacing is equal to A/2. Figure 5.3b and Fig¬ 
ure 5.3c compares the patterns obtained from the T-matrix method versus the point sources. 
The T-matrix beam pattern shows the severe interaction at end fire. The end fire lobe has 3 
dB greater amplitude than the broadside lobe (main lobe). The degradation of amplitudes was 
observed at —132'^, —47^^, 47^^, and 132'^, but the point sources beam pattern goes to adjacent 
nulls at the same angles (Figure 5.3c). 

It was suspected that a linear quadrupolar resonance was responsible that we had not previously 
identified. Exa mi nation of the total scattering amplitude revealed very strong linear quadrupolar 
scattering, about forty times as great as monopolar scattering. An eigenfrequency analysis in the 
vicinity of the operating frequency (474 Hz) was conducted using the COMSOL Multiphysics 
finite-element code. A group of five and only five eigenfrequencies were found very close to 474 
Hz, at about 477 Hz. These could represent modes for which n = 2, m = -2 to 2. Mode shapes for 
these solutions were examined, and all appeared to represent quadrupole modes, with four of the 
five appearing to display predominantly motion in a plane, with quadrupolar symmetry, and with 
exactly one appearing to display linear quadrupolar motion. Confirmation that these most likely 
were indeed quadrupolar resonances was obtained by a second, approximate computation of the 
resonance frequencies of the axially-symmetric eigenfrequencies of a thin elastic shell in water 
[26, pp. 281-284], substituting for the Young’s modulus and Poisson’s ratio the piezoelectric 
elastic compliance modulus values and ^ linear quadrupole eigenfrequency value 

of 473 Hz was found. 

Seven-Element Array Configuration 

Figure 5.4 depicts the 3D normalized far-field beam pattern of the linear array with seven 
evenly-spaced, identical, spherical thin-shell transducer elements at “kfl=l” (/=474 Hz). The 
inter-element spacing is equal to A/2. The width of the main lobe is reduced, and the directiv- 
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ity is enhanced in this configuration compared to the three element array. Two conical nodal 
surfaces are seen between broad side and end fire. 

Figure 5.5 is the normalized vertical beam pattern of the seven evenly-spaced, identical, spher¬ 
ical thin-shell transducers. The inter-element spacing is equal to A/2. Figure 5.5b and Fig¬ 
ure 5.5c compare the patterns obtained from the T-matrix method versus the point sources. Near 
the broadside, the behavior of beam patterns is almost the same as that of the point sources. The 
difference is increased near and at the end fire direction. Separation between lobes cannot be 
detected easily. Three minor lobes were smeared by one wide lobe. While six adjacent nulls 
appear in the beam pattern of the point sources, only one adjacent null occurs in the T-matrix 
method. 


Eight-Element Array Configuration 

Figure 5.6 depicts the 3D normalized far-field beam pattern of the linear array with eight evenly- 
spaced, identical, spherical thin-shell transducer elements at “A:a=l” (/=474 Hz). The inter¬ 
element spacing is equal to A/2. Three conical nodal surfaces were observed. 

Figure 5.7 shows explicitly the difference between the two methods for angles near end fire. 
The T-matrix method does not produce a null to the end fire direction, while the point sources 
do. The Figures 5.7b and 5.7c show that the amplitude values of the T-matrix method away 
from broadside are slightly higher than the point sources in the neighborhood of end fire. 


Three-Element Array Configuration With Different Inter-Element Spacing Distance 

Figure 5.3 shows the severe interaction effect at the end fire where the side lobe has a larger 
amplitude than the main lobe. To investigate further the array interaction of the three-element 
configuration, different inter-element spacing distances for the three-element array were exam¬ 
ined. Figure 5.8 illustrates the effect of an inter-element spacing on the three-element array 
when “A:a=l.” The inter-element spacing distances are selected arbitrarily between A/2-A/4. 
Figures 5.8a and 5.8b show that the A/2 and A/2.2 inter-element spacing distances produce 
larger side lobe than main lobe. Figures 5.8c through 5.8i illustrate the shading of the end fire 
is reduced with respect to the decreasing inter-element spacing distance, and it reaches a min¬ 
ima when spacing distance equals A/3. Although the two methods show some difference in 
amplitude, they have the same beam pattern behavior. 
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Array Configuration With Different Numbers of Elements 

To determine the effect of the number of elements in an array, several different configurations 
were examined. Figures 5.9 and 5.10 represent the far-field beam patterns for differing N (num¬ 
bers of elements in the vertical array). The elements are spaced equally and the inter-element 
spacing is fixed at X/2. For convenience, the configurations with an even number of elements 
and those with an odd number are analyzed in separate figures. 

Figure 5.9 shows the end fire shading is reduced respectively for an increasing number of ele¬ 
ments in an array. When the number of elements in an array is more than three, the nulls which 
are predicted near the end fire in the point sources are shadowed. An analogous shadowing 
effect occurs near the broadside when the number of the elements in an array is equal to 11 
(Figure 5.9e). 

Figure 5.10 shows the configurations for even-numbered array elements. The adjacent null 
is produced at end fire for all even number configuration in point sources due to the axis of 
symmetry while the T-matrix method beam patterns display the array interaction effect, and the 
beam pattern does not terminate at the adjacent nulls at the end fire. 

Figures 5.11 and 5.12 demonstrate the 3D normalized far-field beam pattern for the different 
number configurations corresponding to Figures 5.9 and 5.10. The directionality of the array 
increases with respect to N. The worst directionality appears in Figure 5.11a where the major 
lobe seems to be at the end fire. As shown in Figure 5.11c, Figure 5.Ilf, and the Figure 5.12b 
some of the estimated conical nodal surfaces disappear. 

Summary of Results for Linear Arrays of Spherical Thin-Shell Transducers 

1. The observation of the three-element array configuration shows that for the T-matrix re¬ 
sult, the end-fire lobe is 3 dB greater than the broadside lobe (5.3). The strong array 
interaction was found to be due to strong linear quadrupole scattering. This was traced to 
an uncontrolled (not able to be driven) linear quadrupole resonance of the transducer, at 
477 Hz. The array interaction cause a “shading” in this case and reduce the amplitude of 
the major lobe. 

2. For the T-matrix results, the consistent “shadowing” effect is manifested by a reduction 
in destructive interference, and so true nulls do not occur. 

3. The “shadowing” effect disappears when the inter-element spacing distance is A/3 for 
the three element array (Figure 5.8f). 

4. The configurations of the linear arrays with different N number of the elements, are shown 
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in Figures 5.9 and 5.10, which indicate that the “shading” effect is reduced with respect 
to the increasing number of elements while the “shadowing” effect is increased. It can be 
concluded that the "shading" effect due to scattering is reduced by increasing the number 
of elements in an array. 

5. It should be noted that the beam pattern of both methods is the same at broadside when 
the number of elements in an array is more than three (Figures 5.9 and 5.10). 

6. The least desirable array configuration of the spherical thin-shells was observed for a 
three-element array with A/2 inter element spacing, because the side lobe dominated the 
broadside major lobe (Figure 5.8a). 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=3 and d=> /2 


150 



Figure 5.2: 3D front-left view of normalized far-field beam pattern of vertical linear array with 
three identical spherical thin-shell transducer elements obtained with the T-matrix method,inter¬ 
element spacing d=X/2, /=474 Hz and ka=\. 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2 



(a) (b) (c) 

Figure 5.3: Normalized far-field vertical beam pattern for the modeled linear array with three 
identical spherical thin-shell transducer elements lying along polar axis (6 = 0), inter-element 
spacing d=Xl2, /=474 Hz and ka=\. (a) Orientation of the spherical thin-shell transducers (b) 
A polar plot of the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the 
normalized magnitudes of the vertical far-field beam pattern 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=7 and d=> 12 



Figure 5.4: 3D front-left view of normalized far-field beam pattern of vertical linear array with 
seven identical spherical thin-shell transducer elements obtained with the T-matrix method, 
inter-element spacing d=X/2, /=474 Hz and ka=\ 


NORM FF VERT. BEAM PATTERN in dB N=7 d=>e 



(a) (b) (c) 

Figure 5.5: Normalized far-field vertical beam pattern for the modeled linear array with seven 
identical spherical thin-shell transducer elements lying along polar axis (0 = 0), inter-element 
spacing d=A/2, /=474 Hz and ka=\. (a) Orientation of the spherical thin-shell transducers (b) 
A polar plot of the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the 
normalized magnitudes of the vertical far-field beam pattern 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=8 and d=> 12 



Figure 5.6: 3D front-left view of normalized far-field beam pattern of vertical linear array 
with eight identical spherical thin-shell transducer elements obtained with the T-matrix method, 
inter-element spacing d=X/2, /=474 Hz and ka=\ 


NORM FF VERT. BEAM PATTERN in dB N=8 d=>/2 



(a) (b) (c) 

Figure 5.7: Normalized far-field vertical beam pattern for the modeled linear array with three 
identical spherical thin-shell transducer elements lying along polar axis (0 = 0), inter-element 
spacing <i=A/2, /=474 Hz and ka=\. (a) Orientation of the spherical thin-shell transducers (b) 
A polar plot of the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the 
normalized magnitudes of the vertical far-field beam pattern 
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NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.4 



NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.e 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.8 




(d) d=Xl2.6 


(e) d=X/2.% 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3 



(f) d=X/2 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3.2 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3.4 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/4 





(g) d=X /3.2 (h) d=X /3.4 (i) d=X /4 

Figure 5.8: The comparison of the normalized far-field vertical beam patterns of a modeled 
evenly spaced three-element linear array of spherical thin-shell transducer at /=474 Hz and 
ka=\ with different inter-element spacing distance (a-i) Comparison of polar plot of the nor¬ 
malized vertical far-field beam pattern of T-matrix method against polar plot of the normalized 
vertical far-field beam pattern of the point sources. 
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NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=3 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=5 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=7 and d=>/2 



(a) N=3 


(b) N=5 


(c) N=1 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=9 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=11 and d=>/2 




NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=13 and d=>/2 



I (degree) 


(d) N=9 


(e)N=ll 


(f) iV=13 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=15 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=17 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=19 and d=)/2 




(g) N=15 (h) N=n (i) N=19 


Figure 5.9: The vertical far-field beam patterns for a modeled vertical linear array composed 
of an odd number N of evenly spaced spherical thin-shell transducer elements lying along the 
polar axis (0 = 0) at A/2, /=474 Hz and ka=\ (a-i) Comparison of cartesian plot of the normal¬ 
ized vertical far-field beam pattern of T-matrix method against cartesian plot of the normalized 
vertical far-field beam pattern of the point sources. 
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NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=4 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=6 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=8 and d=>/2 



(a) N=A 


(b) N=6 


(c) N=% 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=10 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=12 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=14 and d=>/2 



(d) N=10 


(e) N=12 


(f) N=U 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=16 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=18 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=20 and d=>/2 



(g) ^=16 


(h)A^=18 


(i) N=20 


Figure 5.10: The vertical far-field beam patterns for a modeled vertical linear array composed 
of an even number N of evenly spaced spherical thin-shell transducer elements lying along the 
polar axis (0 = 0) at A/2, /=474 Hz and ka=\ (a-i) Comparison of cartesian plot of the normal¬ 
ized vertical far-field beam pattern of T-matrix method against cartesian plot of the normalized 
vertical far-field beam pattern of the point sources. 
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Figure 5.11: The 3D normalized far-field beam patterns for a modeled vertical linear array 
composed of an odd number N of evenly spaced spherical thin-shell transducer elements lying 
along the polar axis (0 = 0) at A/2, f-414 Hz and ka=\ (a-i) 3D front-left view of far-field 
beam pattern of vertical linear array 
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Figure 5.12: The 3D normalized far-field beam patterns for a modeled vertical linear array 
composed of an even number N of evenly spaced spherical thin-shell transducer elements lying 
along the polar axis (0 = 0) at A/2, f-414 Hz and ka=\ (a-i) 3D front-left view of far-field 
beam pattern of vertical linear array. 
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5.6.2 The Linear Array with the Ring-Shell Elements 

Three-Element Array Configuration 

Figure 5.13 depicts the 3D normalized far-field beam pattern of a linear array with three evenly- 
spaced, identical, ring-shell transducer elements at “A:a=l” (/=605 Hz). The inter-element spac¬ 
ing is equal to A/2. It is remarkable that the conical nodal surfaces do not degrade relative to 
broadside at end fire 

Figure 5.14 is the normalized vertical beam pattern of the three evenly-spaced, identical, ring- 
shell transducers. The inter-element spacing is again equal to A/2. Figures 5.14b and 5.14c 
compare the beam patterns obtained from the T-matrix method versus the point source method. 
Unlike for the spherical thin shell transducers, the beam patterns for the ring shell transducers 
are nearly the same for both techniques except true nulls do not occur in T-matrix method. This 
is probably due to two things: 

1. The ring-shell transducers are operated at resonance (controlled), where the spherical 
shells were not. 

2. There is no evidence of a spurious (uncontrolled) resonance at the operating frequency 
for the ring-shell transducers, but there was for the spherical shell transducers. 


Seven-Element Array Configuration 

Figure 5.15 depicts the 3D normalized far-field beam pattern of a linear array with seven evenly- 
spaced, identical, ring-shell transducer elements at “A:a=l” (/=605 Hz). The inter-element spac¬ 
ing is again equal to A/2. The width of the main lobe is reduced, and the directivity is enhanced 
in this configuration compared to the three element array. Three conical nodal surfaces can 
clearly be seen in the figure. 

Figure 5.16 is the normalized beam pattern of the seven evenly-spaced, identical, ring-shell 
transducers. The inter-element spacing is equal to A/2. Figures 5.16b and 5.16c compare the 
patterns obtained from the T-matrix method versus the point source calculation. Although the 
beam patterns of the T-matrix seem almost the same as the beam patterns obtained using point 
sources near broadside. Figure 5.16c shows beam patterns obtained using T-matrices do not 
produce nulls as do the point sources. Also, a 2 dB reduction was observed at end fire for the 
T-matrix compared to point sources. 
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Eight-Element Array Configuration 

Figure 5.17 depicts the 3D normalized far-field beam pattern of the linear array with eight 
evenly-spaced, identical, ring-shell transducer elements at “fca=l” (/=605 Hz). The inter¬ 
element spacing is equal to A/2. Four conical nodal surfaces are clearly seen from the figure. 

Although the difference between the two methods is practically indistinguishable at broadside. 
Figure 5.18 shows an interaction effect near the end fire in the detailed plot of amplitudes, 
shown in Figure 5.18c, demonstrating the T-matrix method does not produce nulls at end fire 
while the point source method does. Also it should be noted the adjacent nulls do not exist for 
the T-matrix method. 

Three-Element Array Configuration With Different Inter-Element Spacing Distance 

Figure 5.19 illustrates the effect of an inter-element spacing on the three-element array with 
equal to one. The inter-element spacing distances are arbitrarily selected between A/2-A/4. 
Figures 5.19a through 5.19i show that decreasing the inter-element distances increases array 
interaction. The largest difference is observed for the configuration of the A/3 (Figure 5.19f). 
The difference between the product theory and the T-matrix method are clearly observed near 
the end fire of this configuration. 

Array Configurations With Different Numbers of Elements 

To determine the effect of the number of elements in an array, different configurations were 
examined. Figures 5.20 and 5.21 represent the far-field beam patterns for different numbers N 
of elements. The elements are spaced equally, and inter-element spacing is fixed at A/2. For 
convenience, even numbers of elements and odd numbers of elements are displayed in separate 
figures. 

Figure 5.21 shows the configurations for the even numbers. At end fire, an adjacent null is 
presented with an even number of elements using the point sources due to the axis of symmetry; 
while the beam pattern obtained using the T-matrix shows array interaction effects and produces 
a beam pattern that does not terminate at adjacent nulls at end fire. 

Figures 5.22 and 5.23 demonstrate the 3D normalized far-field beam pattern for the different 
configurations corresponding to figure Figure 5.20 and Figure 5.21. The directionality of the 
array increases with N. Almost all of the estimated conical nodal surfaces by point sources are 
clearly observed. 
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Summary of Results for Linear Arrays of Ring-Shell Transducers 

1. The figures illustrate the difference between the two methods is not remarkable except 
for the end fire direction. Unlike the linear array of spherical thin-shell transducer, the 
‘shading” effect disappears for this configuration. Only the shadowing effect is produced 
near the end fire, and it is remarkable for even numbered configurations. 

2. The most significant scattering effect was observed when the inter-element spacing was 
A/3. A spacing smaller than A/3 resulted in local maxima on axis where the spherical 
shell and point source arrays both show local minima. The relative consistency of axial 
to broadside amplitude for the ring-shell array is probably due to its being generated at 
resonance. 

3. The Figures 5.8 and 5.19 showed that for the spherical thin-shell array, the qualitative 
shape of the beam patterns is the same, e.g., the angles of local maxima and minima. 
However, for the ring-shell array, regardless of the spacing, the on-axis response appears 
to be consistently about 7 dB smaller than the broadside response. 

4. The least desirable array configuration of the ring-shells was observed for A/3 inter 
element-spacing distance (Figure 5.19f). 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=3 and d =>/2 



Figure 5.13: 3D front-left view of normalized far field beam pattern of vertical linear array with 
three ring-shell transducer elements obtained with the T-matrix method,inter-element spacing 
d=X /2, /=605 Hz and ka=\. 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>y2 



(a) (b) (c) 


Figure 5.14: Normalized far-field vertical beam pattern for the modeled linear array with three 
identical ring-shell transducer elements lying along polar axis (6 = 0), inter-element spacing 
d=)i/2, /=605 Hz and ka=\. (a) Orientation of the ring-shell transducers (b) A polar plot of 
the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the normalized 
magnitudes of the vertical far-field beam pattern 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=7 and d =112 



Figure 5.15: 3D front-left view of normalized far field beam pattern of vertical linear array with 
seven ring-shell transducer elements obtained with the T-matrix method, inter-element spacing 
d=X /2, /=605 Hz and ka=\. 


NORM FF VERT. BEAM PATTERN in dB N=7 d=>e 



(a) (b) (c) 


Figure 5.16: Normalized far-field vertical beam pattern for the modeled linear array with seven 
identical ring-shell transducer elements lying along polar axis (0 = 0), inter-element spacing 
d=X/2, /=605 Hz and ka=\. (a) Orientation of the ring-shell transducers (b) A polar plot of 
the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the normalized 
magnitudes of the vertical far-field beam pattern 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=8 and d =>/2 



Figure 5.17: 3D front-left view of normalized far field beam pattern of vertical linear array with 
eight ring-shell transducer elements obtained with the T-matrix method, inter-element spacing 
d=X /2, /=605 Hz and ka=l. 


NORM FF VERT. BEAM PATTERN in dB N=8 d=>/2 



(a) (b) (c) 


Figure 5.18: Normalized far-field vertical beam pattern for the modeled linear array with eight 
identical ring-shell transducer elements lying along polar axis (0 = 0), inter-element spacing 
d=X/2, /=605 Hz and ka=l. (a) Orientation of the ring-shell transducers (b) A polar plot of 
the normalized magnitudes of the vertical far-field beam pattern (c) A plot of the normalized 
magnitudes of the vertical far-field beam pattern 
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NORM FF VERT. BEAM PATTERN in dB N=3 d=>^ NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.2 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.4 





{a)d=Xl2 


(b) d=Xll.l 


(c) d=XllA 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.e NORM FF VERT. BEAM PATTERN in dB N=3 d=>/2.8 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3 





(d) d=Xl2.6 


(e) d=X/2.% 


(f) d=X/2 


NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3.2 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/3.4 NORM FF VERT. BEAM PATTERN in dB N=3 d=>/4 





(g) d=X /3.2 (h) d=X /3.4 (i) d=X /4 

Figure 5.19: The comparison of the far-held vertical beam patterns of a modeled evenly spaced 
three-element linear array of the ring-shell at /=605 Hz and ka=\ with different inter-element 
spacing distance (a-i) Comparison of polar plot of the normalized vertical far-held beam pattern 
of T-matrix method against polar plot of the normalized vertical far-held beam pattern of the 
point sources. 
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NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=3 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=5 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=7 and d=>/2 



(a) N=3 


(b) N=5 


(c) N=1 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=9 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=11 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=13 and d=>/2 



(d) N=9 


(e) N=n 


(f) iV=13 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=15 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=17 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=19 and d=>/2 





(g)N=l5 


(h) N=n 


(i) N=19 


Figure 5.20: The vertical far-field beam patterns for a modeled vertical linear array composed of 
an odd number N of evenly spaced flextensional transducer elements lying along the polar axis 
(0 = 0) at A/2, /=605 Hz and ka-l (a-i) Comparison of cartesian plot of the normalized vertical 
far-field beam pattern of T-matrix method against cartesian plot of the normalized vertical far- 
field beam pattern of the point sources. 
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NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=4 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=6 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=8 and d=>/2 



(a) N=A 


(b) N=6 


(c) N=% 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=10 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=12 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=14 and d=>/2 



(d) N=10 


(e) N=12 


(f) N=U 


NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=16 and d=>/2 NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=18 and d=>/2 




NORM FF VERT. BEAM PATTERN in dB on Z AXiS N=20 and d=>/2 



(g) ^=16 


(h)A^=18 


(i) N=20 


Figure 5.21: The vertical far-field beam patterns for a modeled vertical linear array composed 
of an even number N of evenly spaced spherical thin-shell transducer elements lying along the 
polar axis (0 = 0) at A/2, /=605 Hz and ka=\ (a-i) Comparison of cartesian plot of the normal¬ 
ized vertical far-field beam pattern of T-matrix method against cartesian plot of the normalized 
vertical far-field beam pattern of the point sources. 
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3D Far Field Beam Pattern of Linear Array,ka=1, N=3 and d =>/2 


3D Far Field Beam Pattern of Linear Array,ka=1, N=5 and 


3D Far Field Beam Pattern of Linear Array,ka=1, N=7 and 



d=>/2 



d=>/2 



(b) N=5 


(c) N=1 




Figure 5.22: The 3D normalized far-field beam patterns for a modeled vertical linear array 
composed of an odd number N of evenly spaced spherical thin-shell transducer elements lying 
along the polar axis (0 = 0) at A/2, /=605 Hz and ka-l (a-i) 3D front-left view of normalized 
far-field beam pattern of vertical linear array 
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(h)A^=18 

(g) A^=16 qualitative 



(i) iV=20 


Figure 5.23: The 3D normalized far-field beam patterns for a modeled vertical linear array 
composed of an even number N of evenly spaced spherical thin-shell transducer elements lying 
along the polar axis (0 = 0) at A/2, /=605 Hz and ka=l (a-i) 3D front-left view of normalized 
far-field beam pattern of vertical linear array. 
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5.7 General Conclusion 


The characteristics of beam patterns were obtained using the T-matrix method for modeling a 
vertical linear array comprised of spherical thin shell transducer elements and a vertical linear 
array comprised of ring-shell transducer elements at “fca=l.” Results obtained with the T-matrix 
were compared with for point source arrays. While, both approaches are in agreement near 
broadside, they differ from each other near and at the end fire directions. This is expected 
because the point sources ignore array element interaction while the T-matrix method considers 
this affect. Such interactions, if they stimulate uncontrolled resonances, as is the case for the 
spherical shell transducer, can result in large deviations from point-source behavior. 

Array element interactions cause the degradation or the “shadowing” in amplitude of beam 
pattern amplitudes near end fire. The shadowing or degradation is induced by scattering since 
the scattering produces a destructive interference in an array and changes the radiation away 
from broadside. As a result of these interactions, side lobes smear the minor lobes and shadow 
the adjacent null. In particular, this effect can be seen in 3D beam patterns. The nodal surfaces 
disappear or interfere with the closest side lobes, resulting in the expected nodal surface not 
terminating at an adjacent null (e.g.. Figure 5.4, Figure 5.13). 

To investigate the scattering characteristics, the pressure amplitudes (e.g, outgoing, scattered, 
radiated) are normalized with the monopole radiated pressure amplitude, Fqq. The normalized 
pressure amplitudes are multiplied by phase factors of the far-field Hankel function, in 

order to see the interaction effect on the radiation pattern explicitly (equation (5.5)). The nor¬ 
malized pressure coefficients are listed in Tables 5.3 through 5.8. The scattering characteristics 
can be deduced from these tables. The linear quadrupole component dominates the scattering in 
the array consisting of the spherical shells while monopole scattering dominates the other com¬ 
ponents of the ring-shell array. This is entirely consistent with the previously mentioned facts 
that the ring-shell transducers are operated at (mostly monopole) resonance, where the spher¬ 
ical shell transducers are not, and there exists an uncontrolled linear quadrupole resonance of 
the spherical shell transducers at almost exactly the chosen operating frequency. The differing 
scattering characteristics demonstrate themselves in Figures 5.5 and 5.16. The increasing side 
lobes are observed at end fire of the spherical thin-shells, while the side lobes of the ring-shell 
degrade. 

For further investigation to understand the factors affecting the array interaction, the inter¬ 
element spacing distance was allowed to vary between A/2 —A/4. Resultant figures show 
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the effect of inter-element spacing distances on array sound pressure levels. In particular, when 
the inter-element spacing distance is equal to the ratio of wavelength to the number of elements, 
X/N, significant changes are observed in the beam patterns of Figures 5.19 and 5.8. The ar¬ 
ray interaction can be controlled by adjusting the inter-element spacing distance in order to 
minimize this interaction 


The second factor considered was the affect the number of elements has on beam pattern for 
an array. This array interaction effect is attenuated as the number of elements is increased. In 
conclusion, an optimal number of elements can be determined by analyses of beam patterns to 
enhance or reduce array element interactions. 


N=3 

ka=l f=474 Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.978 

1.000 

0.027 

0.86 

0.00 

146.65 

1 

1,0 

0.068 

0.000 

0.068 

122.39 

90.00 

122.39 

2,0 

1.178 

0.000 

1.178 

320.33 

180.00 

320.33 


0,0 

0.968 

1.000 

0.069 

3.55 

0.00 

119.79 

9 

1,0 

0.000 

0.000 

0.000 

-64.99 

90.00 

-64.99 


2,0 

1.214 

0.000 

1.214 

358.34 

180.00 

358.34 


0,0 

0.978 

1.000 

0.027 

0.86 

0.00 

146.65 

9 

1,0 

0.068 

0.000 

0.068 

-57.61 

90.00 

-57.61 

u 

2,0 

1.178 

0.000 

1.178 

320.33 

180.00 

320.33 


Table 5.3: The normalized pressure coefficient of the linear array consisting of spherical three 
thin-shell transducers 


N=3 

ka=lf=605 Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.855 

1.000 

0.159 

4.05 

0.00 

157.66 

1 

1,0 

0.031 

0.000 

0.031 

10.46 

90.00 

10.46 

2,0 

0.090 

0.104 

0.014 

176.26 

174.59 

343.96 


0,0 

0.918 

1.000 

0.466 

27.73 

0.00 

113.70 

9 

1,0 

0.000 

0.000 

0.000 

-76.86 

90.00 

-76.86 


2,0 

0.090 

0.104 

0.039 

196.55 

174.59 

296.36 


0,0 

0.855 

1.000 

0.159 

4.06 

0.00 

157.60 

9 

1,0 

0.031 

0.000 

0.031 

190.70 

90.00 

190.70 


2,0 

0.090 

0.104 

0.014 

176.27 

174.59 

343.86 


Table 5.4: The normalized pressure coefficient of the linear array consisting of three ring-shell 
transducers 


90 





N=4 

ka=lf=474Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.984 

1.000 

0.031 

1.54 

0.00 

121.45 

1 

1,0 

0.080 

0.000 

0.080 

104.79 

90.00 

104.79 

2,0 

0.885 

0.000 

0.885 

315.54 

180.00 

315.54 


0,0 

0.970 

1.000 

0.054 

2.65 

0.00 

124.43 

9 

1,0 

0.029 

0.000 

0.029 

18.92 

90.00 

18.92 


2,0 

1.249 

0.000 

1.249 

336.80 

180.00 

336.80 


0,0 

0.970 

1.000 

0.054 

2.65 

0.00 

124.43 

9 

1,0 

0.029 

0.000 

0.029 

198.92 

90.00 

198.92 


2,0 

1.249 

0.000 

1.249 

336.80 

180.00 

336.80 


0,0 

0.984 

1.000 

0.031 

1.54 

0.00 

121.45 

4 

1,0 

0.080 

0.000 

0.080 

-75.21 

90.00 

-75.21 


2,0 

0.885 

0.000 

0.885 

315.54 

180.00 

315.54 


Table 5.5: The normalized pressure coefficient of the linear array consisting of four spherical 
thin-shell transducers 


N=4 

ka=lf=605 Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.913 

1.000 

0.196 

10.53 

0.00 

121.42 

1 

1,0 

0.042 

0.000 

0.042 

-17.41 

90.00 

-17.41 

2,0 

0.094 

0.104 

0.016 

182.29 

174.59 

303.89 


0,0 

0.828 

1.000 

0.349 

19.23 

0.00 

128.73 

9 

1,0 

0.032 

0.000 

0.032 

146.43 

90.00 

146.43 


2,0 

0.084 

0.104 

0.030 

187.95 

174.59 

313.81 


0,0 

0.827 

1.000 

0.349 

19.23 

0.00 

128.77 


1,0 

0.032 

0.000 

0.032 

-33.90 

90.00 

-33.90 


2,0 

0.084 

0.104 

0.030 

187.95 

174.59 

313.87 


0,0 

0.914 

1.000 

0.196 

10.54 

0.00 

121.35 

4 

1,0 

0.042 

0.000 

0.042 

162.81 

90.00 

162.81 


2,0 

0.094 

0.104 

0.016 

182.30 

174.59 

303.76 


Table 5.6: The normalized pressure coefficient of the linear array consisting of four ring-shell 
transducers 
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N=8 

ka=l f =474 Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.980 

1.000 

0.032 

1.44 

0.00 

129.79 

1 

1,0 

0.084 

0.000 

0.084 

112.28 

90.00 

112.28 

2,0 

1.033 

0.000 

1.033 

324.51 

180.00 

324.51 


0,0 

0.968 

1.000 

0.058 

2.83 

0.00 

125.04 

9 

1,0 

0.027 

0.000 

0.027 

-22.51 

90.00 

-22.51 


2,0 

1.259 

0.000 

1.259 

345.11 

180.00 

345.11 


0,0 

0.974 

1.000 

0.045 

2.15 

0.00 

126.20 


1,0 

0.033 

0.000 

0.033 

254.49 

90.00 

254.49 


2,0 

1.118 

0.000 

1.118 

329.58 

180.00 

329.58 


0,0 

0.985 

1.000 

0.039 

2.08 

0.00 

113.39 

4 

1,0 

0.051 

0.000 

0.051 

-75.03 

90.00 

-75.03 

‘-r 

2,0 

0.769 

0.000 

0.769 

305.41 

180.00 

305.41 


0,0 

0.985 

1.000 

0.039 

2.08 

0.00 

113.39 


1,0 

0.051 

0.000 

0.051 

104.97 

90.00 

104.97 


2,0 

0.769 

0.000 

0.769 

305.41 

180.00 

305.41 


0,0 

0.974 

1.000 

0.045 

2.15 

0.00 

126.20 


1,0 

0.033 

0.000 

0.033 

74.49 

90.00 

74.49 

VJ 

2,0 

1.118 

0.000 

1.118 

329.58 

180.00 

329.58 


0,0 

0.968 

1.000 

0.058 

2.83 

0.00 

125.04 

7 

1,0 

0.027 

0.000 

0.027 

157.49 

90.00 

157.49 


2,0 

1.259 

0.000 

1.259 

345.11 

180.00 

345.11 


0,0 

0.980 

1.000 

0.032 

1.44 

0.00 

129.79 

8 

1,0 

0.084 

0.000 

0.084 

-67.72 

90.00 

-67.72 

2,0 

1.033 

0.000 

1.033 

324.51 

180.00 

324.51 


Table 5.7: The normalized pressure coefficient of the linear array consisting of eight spherical 
thin-shell transducers 
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N=8 

ka=If=605 Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 


0,0 

0.902 

1.000 

0.178 

8.99 

0.00 

127.67 

1 

1,0 

0.039 

0.000 

0.039 

-13.03 

90.00 

-13.03 

2,0 

0.094 

0.104 

0.015 

180.85 

174.59 

311.18 


0,0 

0.846 

1.000 

0.380 

21.78 

0.00 

124.38 

9 

1,0 

0.024 

0.000 

0.024 

145.83 

90.00 

145.83 


2,0 

0.085 

0.104 

0.032 

190.45 

174.59 

308.74 


0,0 

0.824 

1.000 

0.290 

14.62 

0.00 

134.21 


1,0 

0.020 

0.000 

0.020 

-48.27 

90.00 

-48.27 


2,0 

0.085 

0.104 

0.024 

183.97 

174.59 

320.30 


0,0 

0.841 

1.000 

0.324 

17.69 

0.00 

127.81 

4 

1,0 

0.017 

0.000 

0.017 

125.19 

90.00 

125.19 

‘-r 

2,0 

0.085 

0.104 

0.027 

186.88 

174.59 

312.79 


0,0 

0.841 

1.000 

0.324 

17.69 

0.00 

127.82 


1,0 

0.018 

0.000 

0.018 

-54.94 

90.00 

-54.94 


2,0 

0.085 

0.104 

0.027 

186.88 

174.59 

312.81 


0,0 

0.824 

1.000 

0.290 

14.62 

0.00 

134.19 


1,0 

0.019 

0.000 

0.019 

131.84 

90.00 

131.84 

VJ 

2,0 

0.085 

0.104 

0.024 

183.97 

174.59 

320.27 


0,0 

0.846 

1.000 

0.380 

21.78 

0.00 

124.40 

7 

1,0 

0.025 

0.000 

0.025 

-34.68 

90.00 

-34.68 


2,0 

0.085 

0.104 

0.032 

190.45 

174.59 

308.77 


0,0 

0.902 

1.000 

0.178 

9.00 

0.00 

127.60 

8 

1,0 

0.039 

0.000 

0.039 

167.19 

90.00 

167.19 

2,0 

0.094 

0.104 

0.015 

180.86 

174.59 

311.05 


Table 5.8: The normalized pressure coefficient of the linear array consisting of eight ring-shell 
transducers 
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CHAPTER 6: 

Conclusion and Future Work 


6.1 Conclusion 

This thesis is the first successful practical demonstration of the T-matrix method applied to an 
active sonar array for which a finite-element model was employed to compute the scattering 
properties of a single transducer. 

A piezoelectric thin spherical shell transducer and a class V flextensional transducer were mod¬ 
eled using the COMSOL Multiphysics computer code. Linear array models of the developed 
transducers were created using the T-matrix method. 

Firstly, we modeled a radially polarized piezoelectric spherical thin-shell transducer. Validation 
of the model was done by comparing the computed values of the free-field radiated pressure 
amplitude against those of an approximate analytical model derived from an equivalent circuit 
model of the spherical transducer. The results agreed within 3.7 percent, which is acceptable 
for validation of the finite element code. 

As a second transducer type, a class V flextensional ring-shell transducer was selected in part 
because data are available from a bi-planar billboard array tested in the 1990s [1]. The model 
of the ring-shell flextensional transducer was developed in COMSOL Multiphysics. Computed 
radiation characteristics are in excellent agreement with the manufacturer’s data. The primary 
resonance frequency computed is only 0.8 percent lower than that of the manufacturer’s data, 
and the second resonance computed is only 0.7 percent lower. The calculated source level (SL) 
for a 3000 Vrms drive for the primary resonance frequency was found to be 1.4 dB lower than 
the manufacturer’s data at the primary resonance frequency and 3 dB higher than the manufac¬ 
turer’s data at the secondary resonance frequency. The bandwidth and the mechanical quality 
factor also show excellent agreement with the manufacturer’s data. 

After validation of the model, the scattering T-matrix and radiation coefficients of the developed 
model were computed through quadrupolar order using COMSOL for “A:a” values equal to one. 
This same value was used in the T-matrix method to model the linear array. 

Two linear array models were created. Beam patterns were computed for different configura- 
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tions. The computed beam pattern results were compared against those using point source array 
theory. 

The most significant differences were observed between the two solutions at end fire. This is 
expected because the T-matrix method incorporates the array interaction while the point source 
array theory does not. In particular, we observed that the nulls at end fire are not completed, 
and they are shadowed by scattering. The T- matrix method accounts for the inter-element array 
interaction and simulates the array performance more realistically than the point source array 
theory. 

The T-matrix method can be used as a versatile tool in the design phase of an array with high- 
power active transducer elements. That is because the multi-interaction is accounted for in the 
methodology. The accuracy of the T-matrix method depends upon the finite-element model of 
the single transducer used to create the array. COMSOL Multiphysics produces a good model 
for a single transducer, thereby computing the scattering and radiation characteristics of arrays 
comprised of these transducers. 

It should be noted, the developed finite element model was built for the common case ka=\. 
For higher ka values refinements should be made in the finite element model. Possible im¬ 
provements were investigated for the ring-shell transducer in Chapter 4.4. When the Perfectly 
Matched Layer (PML) boundary condition was applied, the model performed better at frequen¬ 
cies higher than 1600 Hz in determining the SL. 

6.1.1 Future Work 

As a result of the research done on the T-matrix method application for array modeling, several 
other interesting topics for study have arisen. 

In comparison to an approximate method, it is relatively easy to employ the T-matrix method 
in modeling active sonar array element interaction, and so this method should become the pre¬ 
ferred method for active sonar array performance prediction. However, more practical examples 
are needed to demonstrate its applicability and superiority for an actual densely packed, high 
power active sonar array. 

In order to test the accuracy of the T-matrix method, we are encouraged to develop an array 
model of the bi-planar array composed of 24 Sparton of Canada 600 Hz class flextensional 
ring-shell transducers that were tested in the 1990s and compare the data to our computed 
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results using the T-matrix method. 

As an alternative to computation of the mutual acoustic radiation with Pritchard’s approxima¬ 
tion, the development of a mutual acoustic radiation impedance using a T-matrix description 
can be carried out by future research in order to evaluate a projector’s internal stress. 

The accuracy of PML applied in the 3D finite element model can also be investigated for im¬ 
proving the performance of the 3D model at high frequencies. 
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APPENDIX A: 

Transforming the Elasto-Piezoelectric Property Matrix to 
Model a Circumferentially-Polarized Ring as 
Radially-Polarized 


The elasto-piezoelectric property matrix for a piezoelectric ceramic, taking stress T and electric 
field E as independent, is 


Si ■ 


■'ll 


■^13 

0 

0 

0 

0 

0 

d3i 


■ ' 

S2 


A12 

■'ll 

*13 

0 
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Ti 
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■^13 

■^13 

*33 

0 

0 

0 

0 

0 
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T3 
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0 

0 

0 

*44 
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0 

0 

^^15 
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Ss 

= 

0 

0 

0 

0 

*44 

0 
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0 
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Ts 

S6 


0 

0 

0 

0 

0 

*66 

0 

0 
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T6 
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0 

0 

0 

0 

^15 

0 

^11 

0 

0 


El 

Di 


0 

0 

0 

d\5 

0 

0 

0 

^11 
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El 

D3 _ 


. d3i 

d3\ 

d33 
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0 
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0 

0 

^33 J 


.Es . 


(A.1) 


In a circumferentially-polarized ring for a right-handed coordinate basis set, z = 1, r = 2, and 
^ = 3 (axial,radial,circumferential, respectively). Then, denoting / for longitudinal, 5 for shear. 
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(A.2) 


Now, suppose we want to model this ring as radially polarized, while trying to keep its elastic 
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and dielectric matrix as follows. 


Exchange the last two rows and columns (i.e., 8 and 9) 
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(A.3) 


Change the correspondence of the matrix indices with coordinate directions to identify the new 
polarization direction as the “3” direction ((j) = 2, r = 3). Exchanges rows and columns 2 and 3, 
5 and 6, and 8 and 9 
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(A.4) 


The final form is not expected to depend upon whether the transformed polarization direction 
is chosen, r or z. To prove this, we follow the same procedure as above, now exchanging (j) for 
the z polarization direction. We exchange rows and columns (7 and 9). 
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(A.5) 


We change the correspondence of matrix indices with coordinate directions to identify the new 
polarization direction as the “3” direction (we now assigning ^ = I, z = 3) 


S 3 

= Siz 


r 

‘'ll 

‘'12 

‘'13 

0 

0 

0 

^731 

0 

0 


■ 73 

= Tiz ' 

Si 

= Sir 


‘'12 

‘'ll 

'13 

0 

0 

0 

^731 

0 

0 


72 

= Tir 

Si 

= 5/0 


‘'13 

‘'13 

‘'33 

0 

0 

0 

d33 

0 

0 


Ti 

= Ti(i, 

54 

~ Ssz 


0 

0 

0 

S 44 

0 

0 

0 

^7i5 

0 


76 

— Esz 

55 

— Ssr 

= 

0 

0 

0 

0 

‘'44 

0 

0 

0 

^7i5 


T 5 

— Tsr 

56 

~ Ss(j> 


0 

0 

0 

0 

0 

■'66 

0 

0 

0 


74 

= Tsij) 

D 3 

= E>z 


<731 

<731 

d33 

0 

0 

0 

^33 

0 

0 


E 3 

= Ez 

Di 

= Dr 


0 

0 

0 

^7i5 

0 

0 

0 

^11 

0 


El 

= Er 

Di 

= E>,j, 


0 

0 

0 

0 

<7 i 5 

0 

0 

0 

^11 J 


. El 

= E(I> . 


(A.6) 


We exchange rows and columns 1 and 3, 4 and 6, 7 and 9 to put them in regular order: 
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(A.7) 


An alternative procedure to transform the electrical coupling is as follows. Instead of exchang¬ 
ing rows and columns, permute them. Permuting the polarization direction to “r” (i.e., the “2” 
direction): 
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(A.8) 


Now, reidentify coordinate directions and matrix subscripts: 
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Now, put the matrix in regular order: 
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Following the permutation procedure, now make “z” the polarization direction: 


(A.9) 


(A. 10) 
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(A.11) 


*^1 

= Siz 


■'ll 

„E 

■'12 


0 

0 

0 

6?31 

0 

0 


■ n 

= Tiz ■ 


= Sir 


„E 

■'12 

„E 

■'ll 

„E 

^13 

0 

0 

0 

6?31 

0 

0 


Ti 

= Tir 

S 3 

= S/0 


■'IS 

■'IS 

■'ss 

0 

0 

0 

^33 

0 

0 


T 3 

= Ti(j, 

S 4 

~ Ssz 


0 

0 

0 

S 44 

0 

0 

0 

0 

<^15 


T 4 

— Esz 

S 5 

— Ssr 

— 

0 

0 

0 

0 

^44 

0 

0 

^^15 

0 


T 5 

— Tgr 

S6 

II 

-0- 


0 

0 

0 

0 

0 

*66 

0 

0 

0 


Te 

II 

Di 

= Dz 


^31 

^31 

^33 

0 

0 

0 

^33 

0 

0 


El 

= Ez 

Di 

= Dr 


0 

0 

0 

0 

<^15 

0 

0 

^11 

0 


El 

= Er 

Ds 

= D^ 


0 

0 

0 

di5 

0 

0 

0 

0 

F^ 

^11 J 


_E3 

= E^ 


Now, reidentify coordinate directions and matrix subscripts: 
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Now put in regular order: 
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(A.13) 
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APPENDIX B: 
Main Code 


%case: arbitrary number of element on the Z(polar) axis 
clear all 
clc 

element=[2]; 

distance=[2]; 

for numl=l:length(element) 

for num2=l:length(distance) 

clearvars -except element distance numl num2 

clc 

mode=2; % total number of mode 
K2=(mode-l-l)' 2; % monopole,dipole, quadrupole 
% Parameter are adjusted for ka=l 
radius=0.4; %radius of ring-shell [m] 
c=1490; %speed of sound [m/s] 
f=605; 

lambda=c/f; %wavelength [m] 
k=2*pi/lambda; %wave number 
ka=k*radius; % ka=1.0 

d=lambda/distance(num2); %distance between the projector 
kd=k*d; 

%write transducer location for line array 

N=element(numl); %number of projector 

%check the number of element whether odd or even number 

if mod(N,2)==0 %if it is even the spacing along the z axis would 

p=-N/2:N/2-l; 

z=(p*d)-i-d/2; %transducer located through axes(x,y,z) 
else 

p=-(N-l)/2:l:(N-l)/2; 

z=p*d; %transducer located through axes(x,y,z) 
end 
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%other axis would be zero for linear array 

x=zeros(l,N); 

y=zeros(l,N); 

position=[x’ y’ z’]; 

%calculate the relative angle for eaeh transducer 

for varl=l:N 

for var2=l:N 

delta_x=0; 

delta_y=0; 

delta_z=z( var 1 )-z( var2); 

r=sqrt(delta_x. "2+delta_y. ~2+delta_z. ~2); 

theta(varl,var2)=acos(delta_z/r); 

phi(varl,var2)=atan2(delta_x,delta_y); 

end; 

end; 

%Set the diagonal angles zeros(i.e thetall,theta22,theta33...) 
theta(l:N+l:N*N)=0; 

%% Computation of translation Matrix Gij for transducer elements 

% the modes of radiated transducer called as r,m and the modes of the transducer 

% where the radiated pressure defined as incoming pressure reference with 

% it’s origin called as s,u. 

alfa=zeros(K2*N); 

for sph_m=l:N %main ring-shell 

for sph_r=l:N %relative ring-shell 

kr 12=abs(sph_m-sph_r)*kd; 

count2=0; 

for r=0:mode %2nd ring-shell modes r,m 
for m=-r:r 

count2=l-i-count2; %index of Gij matrix column 

countl=0; %index of Gij matrix row 

for s=0:mode %lth ring-shell modes s,u 

for u=-s:s 

count l=count 1-1-1; 

%define the lower limit of summation 
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tO=abs(r-s); 

tl=abs(m-u); 

iftl>tO 

ifmod((tl-t0),2)==0 

tO=tl; 

else 

tO=t 1+1 j 

end; 

end; 

count3=l; 
clear Gtemp3 
for t=t0:2:s+r 

%Van Buren’s paper used for defining the a 
%coeffcient and named this function as "a". 
Gtempl=a(s,t,r,u,m)*shank2(t,krl2); 

%spherical harmonic function calculate sph. 

%harmonic 

M=m-u; %order of Spherical Harmonic 

Gtemp2=sphericalharmonic(t,M,theta(sph_m,sph_r),phi(sph_m,sph_r)); 
%Multiplication of defined spherical harmonic,a 
%coeff and sph. hankel. 

Gtemp3( 1 ,count3)=Gtempl *Gtemp2; 

%increase the count3 t create a Btemplc 
%vector. 

count3=count3+l; 
end; 

%the translation matrix which showed the outgoing wavw 

%from 2 sphere as a incominwave to 1 sphere. 

ifkrl2==0 

G(count 1 ,count2)=0; 

else 

G(countl ,count2)=sum(Gtemp3); 

end; 

end; 
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end; 

countl=l; 

end; 

end; 

%locate the product matrix(T-matrix and Gij) as a submatrix into Alfo matrix 

load(’Tmatrix_ring.mat’) %T-MAtrix computed by the help of COMSOL 

pos 1 a=( 1 +length(G) *(sph_r-1)); 

pos lb=(length(G)+length(G)*(sph_r-1)); 

pos2a= 1 +length(G) * (sph_m-1); 

pos2b=(length(G)+length(G)*(sph_m-l)); 

if sph_r-sph_m==0 

alfa(pos2a:pos2b,posla:poslb)=zeros(length(Tmatrix)); 

else 

alfa(pos2a:pos2b,posla:poslb)=Tmatrix*G; 

end; 

end; 

end; 

%% Calculate Beam Pattern 

clearvars -except alfa N position element distance lambda d k mode numl num2 

load(’ Radiated Pressure Amplitude_ringlV.mat’) %Radiated pressure computed by the help of 

COMSOL 

Pr=repmat(P_radiated,N,l); %repeat the vector for the specified mode for each projector 

%% DEFINE SCATTER VECTOR 

%express in matrix form A*P_s=B 

%A : I-Alfa and B: Free Radiated Pressure multiply alfa 

A=eye(length(alfa))-alfa; 

B=alfa*Pr; 

%Solution for P_s 
P_s=A\B; 

%% Total Pressure Computation 

%Point position,spherical coordinate system according to the top transducer 
sample=30; 

phi=linspace(0,pi,sample); % degree 
%convert radian 
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theta=linspace(0,pi,sample); 

%define the position vectors and unit vectors 
position=position’; 

%Summation 

for ind2=l :length( theta) 

rhat= [cos(phi). * sin(theta(ind2)) ;sin(theta(ind2)). * sin(phi) ;cos(theta(ind2)). *ones( 1 ,length(phi))]; 
for indl=l:length(phi) 

Ptotal=0; 

count=0; 

for proj=l:N %number of projector 

Ptempb=0; 

for n=0:mode 

Ptempa=0; 

for m=-n:n 

count=count+l; 

%Spherical Harmonic Function changed,because theta is not scalar 
%anymore 

Ptempa=Ptempa+(Pr(count)+P_s(count))*sphericalharmonic(n,m,theta(ind2),phi(indl)); % to¬ 
tal Pressure 
end 

Ptempb=Ptempb-i-Ptempa*(li)~(n-i-l); % total Pressure 
end 

Ptotal=Ptotal-i-Ptempb. *exp( 1 i. *k. *(position( 1 ,proj). *rhat( 1 ,ind 1)... 

-l-position(2,proj). *rhat(2,ind 1 )-l-position(3 ,proj). *rhat(3 ,ind 1))); 
end 

Ptotall (indl ,ind2)=Ptotal; 
end; 

%Find the beam pattern 
end; 

%normalization 

H=Ptotall/abs(max(max(Ptotall))); % the max function takes the magnitude of the complex 
value 

b=20*logl0(abs(H)); %Beam Pattern in dB 
shift=50; 
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b =b+shift; 
b(find(b<0)) = 0; 

%% Plots of Results 

X = (b.*(cos(phi)’*sin(theta))); 
y = (b.*(sin(phi)’*sin(theta))); 

z = b.*repmat(cos(theta),sample,!); 

XI =-(b. * (cos(phi) ’ *sin(theta))); 
y 1 =- (b. * (sin(phi) ’ * sin(theta))); 

z 1 =b. *repmat(cos(theta), sample, 1); 

%For convenience in linear Plot convert the angle unit from radian to degree 
phi_d=phi* 180/pi; 
figure l=figure 
axis equal 

set(gca, ’ Fonts ize ’,12) 

surf(x,y,z,b) 

hold on 

surf(xl,yl,zl,b) 

axis square; camlight headlight; 

set(gca, ’ box ’ , ’ off’ ); 

set(gca, ’ xtick’ , [], ’ y tick’ , [], ’ ztick’ , []); 

xlabel(”) 

ylabel(”) 

zlabel(’ ARRAY AXIS’,’FontSize’ , 18) 

title([’3D Far Field Beam Pattern of Linear Array,ka=l, N=’,num2str(N)... 

,’ and d = \lambda/’,num2str(distance(num2))],’FontSize’,15) 

screen=get(0, ’ ScreenSize ’ ); 

end; 

end; 


B.l Code of the Used Functions 

The following code was used in the main code they are written by Joseph L. Day [13] 


116 


shank2 


function hn=shank2(n,x) 

% Computes the Spherical Hankel function of the second kind 
% Written by Joseph L. Day 

hn=sqrt(pi/(2*x))*besselh(n+.5,2,x); 


B.1.1 jprime 

function jp=jprime(n,x) 

% Computes the first derivative of the Spherical Bessel function. 
% Written by Joseph L. Day 

jp=sqrt(pi/(2*x))*(besselj(n-.5,x)-(n+l)*besselj(n+.5,x)/(x)); 


B.1.2 fact 

function y=fact(n) 

% Computes the factorial of n 
% Written by Joseph L. Day 

y=prod(l:n); 


B.1.3 be 

function y=bc(n,m) 

% Computes the binomial coefficient 
% Written by Joseph L. Day 
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ifm<0 

y=0; 

el seif n-m<0 

y=0; 

else 

y=prod(l :n)/(prod(l :m)*prod(l :n-m)); 
end 


B.1.4 a 

function A=a(s,t,r,u,m) 

% computes the function a(s,t,r,u,m) of the Addition Theorem 
% as written in King and Van Buren 

num 1=(2*s+1) * (2*t+1) *fact(s-u) *fact(t-m+u)*fact(r+m) *fact((s+t+r)/2); 
dem 1 =fact((r+t-s)/2) *fact((r+s-t)/2)* fact((s+t-r)/2)*fact(s+t+r+1); 
term 1 =num 1 / dem 1; 

wmin=max([r-s-t,s-r-t-2*m+2*u,t-s-r+2*u])/2; 

wmax=min([s+t-r,r+t-s-2*m+2*u,r+s-t+2*u])/2; 

sum=0; 


for w=wmin:wmax 

sum=sum+(-1) ~ w*bc(s+t-r,(s+t-r)/2+w) *... 

bc(t+r-s,(t+r-s)/2+m-u+w)*bc(s+r-t,(s+r-t)/2-u+w); 

end 

A=(i"(s+t-r))*terml*sum; 
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Initial Distribution List 


1. Defense Technical Information Center 
Ft. Belvoir, Virginia 

2. Dudley Knox Library 
Naval Postgraduate School 
Monterey, California 
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